/* -*- mode: c -*- */

/* Purpose: Given input data in the form of a 2D array, perform a
 * discrete Haar wavelet transform.
 *
 * Author: Ilhan Ozgen-Xian, iozgen@lbl.gov
 */

#include "header.h"

/* Main logic */

int main(int argc, char** argv)
{

  printf("\n\n");

  printf("      __  __    _    __  __ __  __  ___ _____ _   _ \n");
  printf("     |  \\/  |  / \\  |  \\/  |  \\/  |/ _ \\_   _| | | |\n");
  printf("     | |\\/| | / _ \\ | |\\/| | |\\/| | | | || | | |_| |\n");
  printf("     | |  | |/ ___ \\| |  | | |  | | |_| || | |  _  |\n");
  printf("     |_|  |_/_/   \\_\\_|  |_|_|  |_|\\___/ |_| |_| |_|\n\n");

  /* printf("	  +----------------------------+\n"); */
  /* printf("	  |                            |\n"); */
  /* printf("	  |        . ...........       |\n"); */
  /* printf("	  |      ...................   |\n"); */
  /* printf("	  |    .... .. ........... ..  |\n"); */
  /* printf("	  |   ... . ..............  .  |\n"); */
  /* printf("	  |  ..     ..............     |\n"); */
  /* printf("	  |  ..      ..         ..     |\n"); */
  /* printf("	  |          ..         ..     |\n"); */
  /* printf("	  |                            |\n"); */
  /* printf("	  +----------------------------+\n\n"); */

  printf("  I. Ozgen-Xian, G. Kesserwani, D. Caviedes-Voullieme,\n");
  printf("  S. Molins, Z. Xu, D. Dwivedi, J. D. Moulton, C. I. Steefel,\n");
  printf("  Wavelet-based local mesh refinement for rainfall-runoff\n");
  printf("  simulations, Journal of Hydroinformatics, 22, 1059-1077, 2020.\n");
  printf("  doi: 10.2166/hydro.2020.198\n\n");

  

  if (argv[1][0] == '-') {
    
    if (argv[1][1] == 'a') {
      printf("[OK] Carry out multiresolution analysis\n");
      multiresolution_analysis(argc, argv);
    } else if (argv[1][1] == 'm') {
      printf("[OK] Generate mesh\n");
      multiresolution_meshing(argc, argv);
    }
    
  } else {
    fprintf(stderr, "[ERR] Flag not understood: %s\n", argv[1]);
    fprintf(stderr, "[ERR] Usage: ./mammoth -[a|m] <input/file>\n");
  }
  
  return 0;
  
}

int multiresolution_meshing(int argc, char** argv)
{

  FILE *fp;
  FILE *lp;
  FILE *bp;

  char datafilename[255];
  int nlevels;
  double threshold;

  int ncol;
  int nrow;

  double xllcorner;
  double yllcorner;
  double cellsize;
  double nodata;

  int n_bounds;
  int *levels;
  double *boundary_coords;
  double *elevations;
  
  if (argc < 2) {
    fprintf(stderr, "[ERR] No input file provided.\n");
    return -1;
  }

  printf("[OK] Reading in: %s\n", argv[2]);

  fp = fopen(argv[2], "r");

  if (fp == NULL) {
    fprintf(stderr, "[ERR] Couldn't read input file: %s\n", argv[2]);
    return -1;
  }

  char key[1024];

  while (fscanf(fp, "%s", key) != EOF) {

    if (strcmp(key, "datafile") == 0) {
      fscanf(fp, "%s", datafilename);
    } else if (strcmp(key, "N") == 0) {
      fscanf(fp, "%d", &nlevels);
    } else if (strcmp(key, "threshold") == 0) {
      fscanf(fp, "%lf", &threshold);
    } else {
      printf("[!!] Key not recognised: %s\n", key);
      printf("[!!] Continue until known key is encountered\n");
    }

  }

  fclose(fp);

  if (get_metadata(datafilename, &nrow, &ncol, &xllcorner, &yllcorner, &cellsize, &nodata)) {
    return -1;
  }

  printf("[OK] Metadata of %s: nrow: %d, ncol: %d, xll: %lf, yll: %lf, dx: %lf, nodata: %lf\n", datafilename, nrow, ncol, xllcorner, yllcorner, cellsize, nodata);

  /*
   * ----------------------------------------------------------------------
   * read in dem data
   * ----------------------------------------------------------------------
   */

  elevations = calloc(nrow * ncol, sizeof(double));
  get_data(datafilename, nrow, ncol, ncol, elevations);
  printf("[OK] DEM data read\n");

  /*
   * ----------------------------------------------------------------------
   * read in boundary
   * ----------------------------------------------------------------------
   */

  printf("[OK] Reading boundary points {\n");

  bp = fopen("output/boundary.asc", "r");

  if (bp == NULL) {
    fprintf(stderr, "[ERR] No boundary file found. Try running multiresolution analysis first.\n");
    return -1;
  }
  
  int index_x;
  int index_y;
  double z;

  n_bounds = 0;
  char _tmp_char;
  
  while (!feof(bp)) {
    _tmp_char = fgetc(bp);
    if (_tmp_char == '\n') {
      n_bounds ++;
    }
  }

  printf("[ok] Found %d points\n", n_bounds);

  rewind(bp);

  boundary_coords = calloc(3 * n_bounds, sizeof(double));

  for (int counter = 0; counter < n_bounds; counter ++) {
    fscanf(bp, "%d %d %lf", &index_x, &index_y, &z);
    boundary_coords[3*counter] = index_x * cellsize + xllcorner;
    boundary_coords[3*counter+1] = (nrow - index_y) * cellsize + yllcorner;
    boundary_coords[3*counter+2] = z;
#if VERBOSE
    printf("[ok] Read point %d/%d\n", counter + 1, n_bounds);
#endif
  }

  fclose(bp);

#if DEBUG_IO
  FILE *boundary_debug_file;
  boundary_debug_file = fopen("output/_tmp_boundary.asc", "w");
  for (int i = 0; i < n_bounds; i ++) {
    fprintf(boundary_debug_file, "%lf %lf %lf \n", boundary_coords[3*i], boundary_coords[3*i+1], boundary_coords[3*i+2]);
  }
  fclose(boundary_debug_file);
#endif

  printf("[OK] } done\n");

  /*
   * ----------------------------------------------------------------------
   * read in level data
   * ----------------------------------------------------------------------
   */

  int powx = (int) ceil(log2(ncol));
  int powy = (int) ceil(log2(nrow));

  int ncol_padded = ipow(2, powx);
  int nrow_padded = ipow(2, powy);

  levels = calloc(ncol * nrow, sizeof(int));


  for (int i = 0; i < nrow; i ++) {
    for (int j = 0; j < ncol; j ++) {
      levels[j + ncol * i] = (int) nodata;
    }
  }

  int *_tmp_levels = calloc(ncol_padded * nrow_padded, sizeof(int));
  
  lp = fopen("output/l.asc", "r");

  if (lp == NULL) {
    fprintf(stderr, "[ERR] No level file found. Try running multiresolution analysis first.\n");
    return -1;
  }

  while(fscanf(lp, "%s", key) != EOF) {
    
    for (int i = 0; i < nrow_padded; i ++) {
      for (int j = 0; j < ncol_padded; j ++) {
	fscanf(lp, "%d", &_tmp_levels[j + ncol_padded * i]);
      }
    }
    
  }

  fclose(lp);

  for (int i = 0; i < nrow; i ++) {
    for (int j = 0; j < ncol; j ++) {
      levels[j + ncol * i] = _tmp_levels[j + ncol_padded * i];
#if DEBUG_IO
      printf(" %d ", levels[j + ncol * i]);
#endif
    }
#if DEBUG_IO
    printf("\n");
#endif
  }

  printf("[OK] Level data read\n");
  
  /*
   * ----------------------------------------------------------------------
   * build the mesh
   * ----------------------------------------------------------------------
   */

  double a0 = cellsize * cellsize * ipow(2, nlevels-1);

  printf("[OK] Building mesh with largest cell area %.2f m2 {\n\n", a0);

  struct triangulateio in;

  in.numberofpoints = n_bounds;
  in.numberofpointattributes = 1;
  in.numberofsegments = n_bounds;
  in.numberofholes = 0;
  in.numberofregions = 0;
  
  in.pointlist = (double *) calloc(in.numberofpoints * 2, sizeof(double));
  in.pointattributelist = (double *) calloc(in.numberofpoints * in.numberofpointattributes, sizeof(double));
  in.pointmarkerlist = (int *) calloc(in.numberofpoints, sizeof(int));
  in.segmentlist = (int *) calloc(in.numberofsegments * 2, sizeof(int));
  in.segmentmarkerlist = (int *) calloc(in.numberofsegments, sizeof(int));

  /* Each of the columns of this array contains the constraint's x and
     y coordinates are at indices [1] and [2], followed by the
     regional attribute at index [3], followed by the maximum area at
     index [4]. So we have size(regionlist,1)==4. Note that each
     regional attribute is used only if you select the 'A' switch, and
     each area constraint is used only if you select the 'a' switch
     (with no number following), but omitting one of these switches
     does not change the memory layout. Input only, although the
     pointer is copied to the output structure for convenience. */
  in.regionlist = (double *) calloc(in.numberofregions * 4, sizeof(double));

  for (int i = 0; i < in.numberofpoints; ++ i) {
    in.pointmarkerlist[i] = 1;
    //in.pointattributelist[i] = boundary_coords[3 * i + 2]; // set z coordinates
    for (int j = 0; j < 2; ++ j) {
      in.pointlist[2 * i + j] = boundary_coords[3 * i + j];
    }
  }

  for (int i = 0; i < in.numberofsegments-1; ++ i) {
    in.segmentlist[2*i] = in.numberofsegments-i;
    in.segmentlist[2*i+1] = in.numberofsegments-1-i;
  }

  in.segmentlist[2*in.numberofpoints-2] = 1;
  in.segmentlist[2*in.numberofpoints-1] = in.numberofsegments;

  struct triangulateio out[nlevels];

  out[0].pointlist = (double *) NULL;
  out[0].pointattributelist = (double *) NULL;
  out[0].pointmarkerlist = (int *) NULL;
  out[0].trianglelist = (int *) NULL;
  out[0].triangleattributelist = (double *) NULL;
  out[0].neighborlist = (int *) NULL;
  out[0].segmentlist = (int *) NULL;
  out[0].segmentmarkerlist = (int *) NULL;

  printf("[ok] Initial triangulation...");
  
  char triangleargs[1024];
  sprintf(triangleargs, "QDpa%f", a0);
  triangulate(triangleargs, &in, &out[0], (struct triangulateio *) NULL);

  printf(" done\n");

  write_avs(&out[0], 1, elevations, ncol, nrow, xllcorner, yllcorner, cellsize, nodata);

  for (int current_level = 1; current_level < nlevels; ++ current_level) {

    out[current_level-1].trianglearealist = (double *) calloc(out[current_level-1].numberoftriangles, sizeof(double));

    // tag the stuff
    for (int i = 0; i < out[current_level-1].numberoftriangles; ++ i) {

      int p1 = out[current_level-1].trianglelist[3*i] - 1;
      int p2 = out[current_level-1].trianglelist[3*i+1] - 1;
      int p3 = out[current_level-1].trianglelist[3*i+2] - 1;

      double x1 = out[current_level-1].pointlist[2*p1] - xllcorner;
      double y1 = out[current_level-1].pointlist[2*p1+1] - yllcorner;
    
      double x2 = out[current_level-1].pointlist[2*p2] - xllcorner;
      double y2 = out[current_level-1].pointlist[2*p2+1] - yllcorner;

      double x3 = out[current_level-1].pointlist[2*p3] - xllcorner;
      double y3 = out[current_level-1].pointlist[2*p3+1] - yllcorner;

      double xmin = min(x1, min(x2, x3));
      double ymin = min(y1, min(y2, y3));
      double xmax = max(x1, max(x2, x3));
      double ymax = max(y1, max(y2, y3));

      int ixmin = (int) floor(xmin / cellsize);
      int iymin = (int) floor(ymin / cellsize);
      int ixmax = (int) floor(xmax / cellsize);
      int iymax = (int) floor(ymax / cellsize);

      /* double xc = (x1 + x2 + x3) / 3.0; */
      /* double yc = (y1 + y2 + y3) / 3.0; */
      /* double a  = 0.5 * (- y2 * x3 + y1 * (-x2 + x3) + x1 * (y2 - y3) + x2 * y3); */
      /* double c  = 1.0 / (2.0 * a); */

      /* double f1 = y1 * x2; */
      /* double f2 = x1 * y3; */
      /* double f3 = y3 - y1; */
      /* double f4 = x1 - x3; */
    
      /* double g1 = x1 * y2; */
      /* double g2 = y1 * x2; */
      /* double g3 = y1 - y2; */
      /* double g4 = x2 - x1; */

      int flag = -1;

      /* printf("[ok] Triangle #%d/%d\n", i+1, mid.numberoftriangles); */
      /* printf("\t %d %d %d / %d\n", p1, p2, p3, mid.numberofpoints); */
      /* printf("\t %d %d : %d %d\n", ixmin, ixmax, iymin, iymax); */
      /* printf("\t %f %f : %f %f\n", xmin, xmax, ymin, ymax); */
      
      for (int iy = iymin; iy < iymax; ++ iy) {
	for (int ix = ixmin; ix < ixmax; ++ ix) {
	  /* printf("%d < %d?\n", flag, levels[ix + ncol * iy]); */
	  flag = (flag < levels[ix + ncol * (nrow - iy)]) ? levels[ix + ncol * (nrow -iy)] : flag;
	}
      }

      /* printf("\t flag: %d...", flag); */

      if (flag >= current_level) {
	/* printf(" tagged\n"); */
	out[current_level-1].trianglearealist[i] = a0 / ipow(2, current_level);
      }

    }
  
    out[current_level].pointlist = (double *) NULL;
    out[current_level].pointattributelist = (double *) NULL;
    out[current_level].pointmarkerlist = (int *) NULL;
    out[current_level].trianglelist = (int *) NULL;
    out[current_level].triangleattributelist = (double *) NULL;
    out[current_level].segmentlist = (int *) NULL;
    out[current_level].segmentmarkerlist = (int *) NULL;

    printf("[ok] Triangulation #%d ...", current_level + 1);

    triangulate("QDarp", &out[current_level-1], &out[current_level], (struct triangulateio *) NULL);
    
    printf(" done\n");

    write_avs(&out[current_level], current_level+1, elevations, ncol, nrow, xllcorner, yllcorner, cellsize, nodata);

  }
  
  printf("[OK] } done\n");


  free(in.pointlist);
  free(in.pointattributelist);
  free(in.pointmarkerlist);
  free(in.regionlist);

  free(out[0].pointlist);
  free(out[0].pointattributelist);
  free(out[0].pointmarkerlist);
  free(out[0].trianglelist);
  free(out[0].triangleattributelist);
  free(out[0].neighborlist);
  
  free(levels);
  free(elevations);
  free(boundary_coords);
  
  return 0;
  
}

int multiresolution_analysis(int argc, char** argv)
{

  FILE *fp;

  char datafilename[255];
  int nlevels;
  double threshold;

  int ncol;
  int nrow;

  double xllcorner;
  double yllcorner;
  double cellsize;
  double nodata;

  double *data;

  if (argc < 2) {
    fprintf(stderr, "[ERR] No input file provided.\n");
    return -1;
  }

  printf("[OK] Reading in: %s\n", argv[2]);

  fp = fopen(argv[2], "r");

  if (fp == NULL) {
    fprintf(stderr, "[ERR] Couldn't read input file: %s\n", argv[2]);
    return -1;
  }

  char key[1024];

  while (fscanf(fp, "%s", key) != EOF) {

    if (strcmp(key, "datafile") == 0) {
      fscanf(fp, "%s", datafilename);
    } else if (strcmp(key, "N") == 0) {
      fscanf(fp, "%d", &nlevels);
    } else if (strcmp(key, "threshold") == 0) {
      fscanf(fp, "%lf", &threshold);
    } else {
      printf("[!!] Key not recognised: %s\n", key);
      printf("[!!] Continue until known key is encountered\n");
    }

  }

  fclose(fp);

  if (get_metadata(datafilename, &nrow, &ncol, &xllcorner, &yllcorner, &cellsize, &nodata)) {
    return -1;
  }

  printf("[OK] Metadata of %s: nrow: %d, ncol: %d, xll: %lf, yll: %lf, dx: %lf, nodata: %lf\n", datafilename, nrow, ncol, xllcorner, yllcorner, cellsize, nodata);

  int powx = (int) ceil(log2(ncol));
  int powy = (int) ceil(log2(nrow));

  int ncol_padded = ipow(2, powx);
  int nrow_padded = ipow(2, powy);

  printf("[OK] Padded data to %d-by-%d\n", nrow_padded, ncol_padded);

  /* read in data */
  data = (double*) calloc(ncol_padded * nrow_padded, sizeof(double));
  for (int i = 0; i < nrow_padded; i ++) {
    for (int j = 0; j < ncol_padded; j ++) {
      data[j + ncol_padded * i] = nodata;
    }
  }

  if (get_data(datafilename, nrow, ncol, ncol_padded, data)) {
    return -1;
  }

  /*
   * ----------------------------------------------------------------------
   * Compression step
   * ----------------------------------------------------------------------
   */

  double *m = (double*) calloc(ncol_padded * nrow_padded, sizeof(double));
  double *_data = (double*) calloc(ncol_padded * nrow_padded, sizeof(double));
  double *_m = (double*) calloc(ncol_padded * nrow_padded, sizeof(double));

  memcpy(_data, data, ncol_padded * nrow_padded * sizeof(double));

  int _nrow = nrow_padded;
  int _ncol = ncol_padded;

  for (int n = 0; n < nlevels; n ++) {

    /*
     * If non-divisible by 2, DHT can't be applied since it depends on
     * pair-wise comparisons. In this case, we exit with fail.
     */
    if (_ncol % 2 != 0 || _nrow % 2 != 0) {
      fprintf(stderr, "[ERR] %d-by-%d matrix can't be compressed by DHT.\n", _ncol, _nrow);
      return -1;
    }

    /*
     * carry out the DHT
     */

    compress(_data, _m, _nrow, _ncol);

    /*
     * Now get the new image matrix. We know that it is located at the
     * upper left corner of the compressed image M, and has dimensions
     * ncol/2 x nrow/2.
     */

    _nrow = _nrow / 2;
    _ncol = _ncol / 2;

    /* Resize _data */
    _data = (double*) realloc(_data, _ncol * _nrow * sizeof(double));

    printf("[OK] Data is resized\n");

    /* Transfer _data */
    for (int i = 0; i < _nrow; i ++) {
      for (int j = 0; j < _ncol; j ++) {
	_data[j + i * _ncol] = _m[j + i * 2 * _ncol];
      }
    }

    printf("[OK] Data is transferred\n");

    /* Transfer m */
    for (int i = 0; i < 2 * _nrow; i ++) {
      for (int j = 0; j < 2 * _ncol; j ++) {
	m[j + i * ncol_padded] = _m[j + i * 2 * _ncol];
      }
    }

#if DEBUG_IO

    char filename[1024];
    sprintf(filename, "output/m-%d.asc", n);

    FILE *fp = fopen(filename, "w");

    for (int i = 0; i < 2 * _nrow; i ++) {
      for (int j = 0; j < 2 * _ncol; j ++) {
	fprintf(fp, "%f ", _m[j + i * 2 * _ncol]);
      }
      fprintf(fp, "\n");
    }

    fclose(fp);

#endif

    printf("[OK] M is transferred\n");

    /* Resize _m */
    _m = (double*) realloc(_m, _ncol * _nrow * sizeof(double));

    printf("[OK] M is resized\n");

  }

  /*
   * ----------------------------------------------------------------------
   * Decompression step
   * ----------------------------------------------------------------------
   * Decompressing the DHT, in each level we look at the details and
   * see if they are significant.  i = 0 is the coarsest level, and
   * the subsequent finer levels are i = 1, 2, ... , N - 1,
   * respectively.
   * ----------------------------------------------------------------------
   */

  int *levels = (int*) calloc(nrow_padded * ncol_padded, sizeof(int));

  /* Set resolution level to highest at the beginning */
  for (int i = 0; i < nrow_padded; i ++) {
    for (int j = 0; j < ncol_padded; j ++) {
      levels[j + ncol_padded * i] = nlevels - 1;
    }
  }

  _nrow = nrow_padded / ipow(2, nlevels - 1);
  _ncol = ncol_padded / ipow(2, nlevels - 1);

  _m = (double*) realloc(_m, _ncol * _nrow * sizeof(double));

  for (int i = 0; i < _nrow; i ++) {
    for (int j = 0; j < _ncol; j ++) {
      _m[j + _ncol * i] = m[j + i * ncol_padded];
    }
  }

  int *flag = (int*) calloc(_nrow * _ncol / 4, sizeof(double));

  for (int n = 0; n < nlevels; n ++) {

    printf("[OK] Thresholding image at level %d\n", n);

    /* Perform hard thresholding to predict refinement levels */
    hard_threshold(_m, flag, n, nlevels, _nrow, _ncol, threshold);

    /* Tag cells that are sufficiently refined */
    for (int i = 0; i < _nrow/2; i ++) {
      for (int j = 0; j < _ncol/2; j ++) {

	if (flag[j + i * _ncol/2] == 1) {

	  int coef = ipow(2, nlevels - n);

	  if (levels[coef * j + coef * i * ncol_padded] > n) {

	    for (int m = coef * i; m < coef * i + coef; m ++) {
	      for (int p = coef * j; p < coef * j + coef; p ++) {
		levels[p + m * ncol_padded] = n;
	      }
	    }

	  }

	}

      }
    }

#if DEBUG_IO

    char filename[1024];
    sprintf(filename, "f-%d.asc", n);

    FILE *fp = fopen(filename, "w");

    for (int i = 0; i < _nrow/2; i ++) {
      for (int j = 0; j < _ncol/2; j ++) {
	fprintf(fp, "%d ", flag[j + i * _ncol/2]);
      }
      fprintf(fp, "\n");
    }

    fclose(fp);

#endif

    /* Decompress the blur image to move to the next refinement level */
    decompress(_m, _nrow, _ncol);

    /* Transfer m */
    for (int i = 0; i < _nrow; i ++) {
      for (int j = 0; j < _ncol; j ++) {
	m[j + i * ncol_padded] = _m[j + i * _ncol];
      }
    }

    printf("[OK] M is transferred\n");

    /*
     * Adjust the dimensions. We are now refining so _ncol and _nrow
     * need to be doubled
     */
    _ncol = 2 * _ncol;
    _nrow = 2 * _nrow;

    /* Resize _m */
    _m = (double*) realloc(_m, _ncol * _nrow * sizeof(double));

    printf("[OK] M is resized\n");

    /* Transfer m back */
    for (int i = 0; i < _nrow; i ++) {
      for (int j = 0; j < _ncol; j ++) {
	_m[j + i * _ncol] = m[j + i * ncol_padded];
      }
    }

    printf("[OK] M is transferred back\n");

    /* Resize flag */
    flag = (int*) realloc(flag, (_nrow * _ncol / 4) * sizeof(double));

    printf("[OK] Flag is resized\n");

  }

  /*
   * ----------------------------------------------------------------------
   * I/O
   * ----------------------------------------------------------------------
   */

#if DEBUG_IO
  FILE *mptr; /* Terrain data */
  mptr = fopen("output/m.asc", "w");
#endif

  FILE *lptr; /* Level data   */
  lptr = fopen("output/l.asc", "w");

  for (int i = 0; i < nrow_padded; i ++) {
    for (int j = 0; j < ncol_padded; j ++) {
#if DEBUG_IO
      fprintf(mptr, "%f ", m[j + ncol_padded * i]);
#endif
      fprintf(lptr, "%d ", levels[j + ncol_padded * i]);
    }
#if DEBUG_IO
    fprintf(mptr, "\n");
#endif
    fprintf(lptr, "\n");
  }

  /*
   * ----------------------------------------------------------------------
   * Meshing
   * ----------------------------------------------------------------------
   */

  /* Resize _data */
  _data = (double*) realloc(_data, ncol * nrow * sizeof(double));

  for (int i = 0; i < nrow; i ++) {
    for (int j = 0; j < ncol; j ++) {
      _data[j + ncol * i] = data[j + ncol_padded * i];
    }
  }

  #if DEBUG_IO
  FILE *_tmp;
  _tmp = fopen("_tmp.asc", "w");

  for (int i = 0; i < nrow; i ++) {
    for (int j = 0; j < ncol; j ++) {
      fprintf(_tmp, "%f ", _data[j + ncol * i]);
    }
    fprintf(_tmp, "\n");
  }
  fclose(_tmp);
  #endif

  printf("[OK] Get boundary using square tracing... ");
  
  get_border(_data, nodata, nrow, ncol);

  printf("done\n");

/*   FILE *bp; */
/*   FILE *op; */

/*   bp = fopen("output/boundary.asc", "r"); */

/*   int index_x; */
/*   int index_y; */
/*   double z; */

/*   int n_bounds = 0; */
/*   char _tmp_char; */
  
/*   while (!feof(bp)) { */
/*     _tmp_char = fgetc(bp); */
/*     if (_tmp_char == '\n') { */
/*       n_bounds ++; */
/*     } */
/*   } */

/*   printf("[ok] Found %d points\n", n_bounds); */

/*   rewind(bp); */

/*   double *boundary_coords = calloc(3 * n_bounds, sizeof(double)); */

/*   for (int counter = 0; counter < n_bounds; counter ++) { */
/*     fscanf(bp, "%d %d %lf", &index_x, &index_y, &z); */
/*     boundary_coords[3*counter] = index_x * cellsize + xllcorner; */
/*     boundary_coords[3*counter+1] = (nrow - index_y) * cellsize + yllcorner; */
/*     boundary_coords[3*counter+2] = z; */
/* #if VERBOSE */
/*     printf("[ok] Read point %d/%d\n", counter + 1, n_bounds); */
/* #endif */
/*   } */

/*   fclose(bp); */

/*   op = fopen("output/boundary.poly", "w"); */

/*   fprintf(op, "%d %d %d %d\n", n_bounds, 2, 1, 0); */

/*   for (int i = 0; i < n_bounds; ++ i) { */
/*     fprintf(op, "%d %f %f %f\n", i+1, boundary_coords[3*i], boundary_coords[3*i+1], boundary_coords[3*i+2]); */
/*   } */

/*   fprintf(op, "%d %d\n", n_bounds, 0); */

/*   for (int i = 0; i < n_bounds-1; ++ i) { */
/*     fprintf(op, "%d %d %d\n", i+1, n_bounds-i, n_bounds-1-i); */
/*   } */

/*   fprintf(op, "%d %d %d\n0 0\n", n_bounds, 1, n_bounds); */

/*   fclose(op); */

/*   free(boundary_coords); */

  //FILE *dp;
  //
  //dp = fopen("output/boundary.asc", "r");
  //
  //int index_x;
  //int index_y;
  //double _tmp;
  //
  //while (!feof(dp)) {
  //  
  //  fscanf(dp, "%d %d %lf", &index_x, &index_y, &_tmp);
  //  struct node p = {.x = index_x * cellsize + xllcorner, .y = index_y * cellsize + yllcorner, .next = NULL};
  //  push(&bcpoints, &p);
  //  
  //}
  //
  //fclose(dp);
  //
  //printf("[OK] Boundary polygon constructed with %ld points\n", bcpoints.len);
  //
  //struct node *tmp = pop(&bcpoints);
  //
  //printf("[OK] %lf %lf\n", tmp->x, tmp->y);
  //
  /* size_t len_xybc = bcpoints.len; */
  /* xybc = (double *) calloc(2 * len_xybc, sizeof(double)); */

  /* printf("[OK] Processing points: "); */
  /* for (size_t i = 0; i < len_xybc; i ++) { */
  /*   printf(" %ld ", i); */
  /*   struct node *p = pop(&bcpoints); */
  /*   xybc[2*i]   = p->x * cellsize + xllcorner; */
  /*   xybc[2*i+1] = p->y * cellsize + yllcorner; */
  /* } */
  /* printf("\n"); */

  /*  */

  /*
   * ----------------------------------------------------------------------
   * Housekeeping
   * ----------------------------------------------------------------------
   */

#if DEBUG_IO
  fclose(mptr);
#endif

  fclose(lptr);

  /* free(xybc); */
  free(levels);
  free(flag);
  free(_data);
  free(m);
  free(data);

  printf("[OK] Finished successfully\n");

  return 0;

}

/*
 * Functions related to reading in the data set
 */

int get_metadata(char *datafilename, int *nrow, int *ncol, double *xllcorner, double *yllcorner, double *cellsize, double *nodata)
{

  FILE *dp;

  printf("[OK] Reading in: %s\n", datafilename);
  dp = fopen(datafilename, "r");

  if (dp == NULL) {
    fprintf(stderr, "[ERR] Couldn't read data file: %s\n", datafilename);
    return -1;
  }

  char key[1024];

  while (fscanf(dp, "%s", key) != EOF) {

    if (strcmp(key, "ncols") == 0) {
      fscanf(dp, "%d", ncol);
    } else if (strcmp(key, "nrows") == 0) {
      fscanf(dp, "%d", nrow);
    } else if (strcmp(key, "xllcorner") == 0) {
      fscanf(dp, "%lf", xllcorner);
    } else if (strcmp(key, "yllcorner") == 0) {
      fscanf(dp, "%lf", yllcorner);
    } else if (strcmp(key, "cellsize") == 0) {
      fscanf(dp, "%lf", cellsize);
    } else if (strcmp(key, "nodata_value") == 0) {
      fscanf(dp, "%lf", nodata);
      return 0;
    } else {
      printf("[!!] Key not recognised: %s\n", key);
      printf("[!!] Continue until known key is encountered\n");
    }

  }

  fclose(dp);

  return 0;
}

int get_data(char *datafilename, int nrow, int ncol, int ncol_padded, double *data)
{

  FILE *dp;
  dp = fopen(datafilename, "r");

  if (dp == NULL) {
    fprintf(stderr, "[ERR] Couldn't read data file: %s\n", datafilename);
    return -1;
  }

  double nodata;
  char key[1024];

  while (fscanf(dp, "%s", key) != EOF) {

    if (strcmp(key, "nodata_value") == 0) {

      fscanf(dp, "%lf", &nodata);
#if VERBOSE
      printf("[OK] Starting to read data\n");
#endif
      
      for (int i = 0; i < nrow; i ++) {
	for (int j = 0; j < ncol; j ++) {
	  fscanf(dp, "%lf", &data[j + ncol_padded * i]);
	}
      }

    }

  }

  return 0;

}

/*
 * Functions related to discrete Haar wavelet transform
 */

int compress(double* data, double *m, int nrow, int ncol)
{
  /*
   *      | B : V |
   * M =  |-------|
   *      | H : D |
   *
   * where B is the blur image, V is the vertical detail, H is the
   * horizontal detail, and D is the diagonal detail.
   */

  printf("[OK] Compressing a %d-by-%d matrix\n", nrow, ncol);

  /*
   * Get the coefficient matrix for left matrix multiplication:
   *
   * C = Wm A
   *
   * This matrix needs to have dimension `nrow x nrow'. After that,
   * carry out the matrix multiplication to get the column-compressed
   * matrix C.
   */

  double *wm;
  wm = (double*) calloc(nrow * nrow, sizeof(double));

  get_wavelet_matrix(wm, nrow);

  double *c;
  c = (double*) calloc(nrow * ncol, sizeof(double));

  matmul(wm, nrow, nrow, data, nrow, ncol, c);

  /*
   * Carry out the row compression. We need a coefficient matrix with
   * dimension `ncol x ncol' to multiply:
   *
   * M = C Wn^T
   *
   * This will compress the image one level.
   */

  double *wn;
  wn = (double*) calloc(ncol * ncol, sizeof(double));

  get_wavelet_matrix(wn, ncol);
  transpose(wn, ncol, ncol);

  matmul(c, nrow, ncol, wn, ncol, ncol, m);

  free(wn);
  free(c);
  free(wm);

  return 0;
}

int decompress(double *data, int nrow, int ncol)
{

  printf("[OK] Decompressing a %d-by-%d matrix\n", nrow, ncol);

  double *wm = (double*) calloc(nrow * nrow, sizeof(double));
  double *wn = (double*) calloc(ncol * ncol, sizeof(double));

  double *c1 = (double*) calloc(nrow * ncol, sizeof(double));
  double *c2 = (double*) calloc(nrow * ncol, sizeof(double));

  /* Reconstruct the DHT to move to finer resolution */

  get_wavelet_matrix(wm, nrow);
  get_wavelet_matrix(wn, ncol);

  transpose(wm, nrow, nrow);

  matmul(wm, nrow, nrow, data, nrow, ncol, c1);
  matmul(c1, nrow, ncol, wn, ncol, ncol, c2);

  for (int i = 0; i < nrow; i ++) {
    for (int j = 0; j < ncol; j ++) {
      data[j + ncol * i] = c2[j + i * ncol];
    }
  }

  free(wm);
  free(wn);
  free(c1);
  free(c2);

  return 0;
}

int hard_threshold(double *m, int *flag, int n, int nlevels, int nrow, int ncol, double threshold)
{
  /*
   * Hard thresholding as explained in (Caviedes-Voullieme and
   * Kesserwani, 2015: Benchmarking a multi-resolution discontinuous
   * Galerkin shallow water model, Advances in Water Resources,
   * 86:14-31).
   */

  /*
   *      | B : V |
   * M =  |-------|
   *      | H : D |
   */

  double _b;
  double _v;
  double _h;
  double _d;

  double _err = ipow(2, n - nlevels) * threshold;

  for (int i = 0; i < nrow/2; i ++) {
    for (int j = 0; j < ncol/2; j++) {

      double max_b = 0.0;
      double max_v = 0.0;
      double max_h = 0.0;
      double max_d = 0.0;

      double max_hvd = 0.0;

      double abs_b = 0.0;
      double abs_v = 0.0;
      double abs_h = 0.0;
      double abs_d = 0.0;

      _b = m[j + ncol * i];
      _v = m[(j + ncol/2) + ncol * i];
      _h = m[j + ncol * (nrow/2 + i)];
      _d = m[(j + ncol/2) + ncol * (nrow/2 + i)];

      /* printf(">>> %d %d %f %f %f %f\n", j, i, _b, _v, _h, _d); */

      abs_b = fabs(_b);
      max_b = (max_b < abs_b) ? abs_b : max_b;

      abs_v = fabs(_v);
      max_v = (max_v < abs_v) ? abs_v : max_v;

      abs_h = fabs(_h);
      max_h = (max_h < abs_h) ? abs_h : max_h;

      abs_d = fabs(_d);
      max_d = (max_d < abs_d) ? abs_d : max_d;

      max_hvd = (max_v > max_h) ? max_v : max_h;
      max_hvd = (max_hvd > max_d) ? max_hvd : max_d;

      /* normalise */
      if (max_b > 1.0) {
	max_hvd = max_hvd / max_b;
      }

      if (max_hvd < _err) {
	/* Can be coarsened */
	flag[j + (ncol/2) * i] = 1;
	/* printf(">>> coarsen: %d %d %f %f %f %f\n", j, i, max_v, ); */
      } else {
	/* Can't be coarsened */
	flag[j + (ncol/2) * i] = -1;
      }

    }
  }

  return 0;
}