/*@@
   @file      ieee_extract.c
   @date      Fri 14 Dec 2001
   @author    Thomas Radke
   @desc
              Utility program to extract a hyperslab from an IEEEIO datafile.
   @enddesc
   @version   $Id$
 @@*/

#define  MAXDIM    3
#define  MAXNAMESIZE  100

#include <stdio.h>
#include <stdlib.h>

/* CCTK includes */
#include "cctk.h"

/* FlexIO includes */
#include "IOProtos.h"
#include "IEEEIO.h"

/* the rcs ID and its dummy function to use it */
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusPUGHIO_IOFlexIO_util_ieee_extract_c)


typedef struct
{
  double x, y, z;
} double3_t;


int main (int argc, char **argv)
{
  int i, j;
  IOFile infile, outfile;
  int nDatasets, nAttributes;
  int datatype;
  int rank;
  int dims [MAXDIM], chunk_dims [MAXDIM], chunk_origin [MAXDIM];
  int attrType, attrLen;
  char attrName [MAXNAMESIZE];
  void *data, *attrData;
  double3_t delta, min_ext, *max_ext;


  if (argc != 4)
  {
    printf ("Usage: %s <region> <infile> <outfile>\n", argv [0]);
    printf ("   eg. %s 128x128x128+32+32+32 foo.ieee bar.ieee\n", argv [0]);
        printf ("       selects a 128-cubed region with origin (32, 32, 32)\n");
    return (0);
  }

  if (sscanf (argv [1], "%dx%dx%d+%d+%d+%d",
                &chunk_dims [0], &chunk_dims [1], &chunk_dims [2],
                &chunk_origin [0], &chunk_origin [1], &chunk_origin [2]) != 6)
  {
    printf ("ERROR: Cannot parse region argument '%s'\n\n", argv [1]);
    return (-1);
  }
  printf ("Selecting region (%d, %d, %d) at origin (%d, %d, %d)\n",
            chunk_dims [0], chunk_dims [1], chunk_dims [2],
            chunk_origin [0], chunk_origin [1], chunk_origin [2]);

  infile = IEEEopen (argv [2], "r");
  if (! IOisValid (infile))
  {
    printf ("ERROR: Could not open input file '%s'\n\n", argv [2]);
    return (1);
  }
  nDatasets = IOnDatasets (infile);
  printf ("Input file '%s' contains %d datasets\n", argv [2], nDatasets);

  outfile = IEEEopen (argv [3], "w");
  if (! IOisValid (outfile))
  {
    printf ("Could not open output file '%s'\n", argv [3]);
    return (1);
  }

  for (i = 0; i < nDatasets; i++)
  {
    if (IOreadInfo (infile, &datatype, &rank, dims, MAXDIM) <= 0)
    {
      printf ("ERROR: Cannot read info of datatset %d, "
              "ignoring dataset...\n\n", i);
      continue;
    }

    nAttributes = IOnAttributes (infile);
    printf ("Processing dataset %d with %d attributes\n", i, nAttributes);
    printf ("Dataset dimensions are (%d, %d, %d)\n",
                dims [0], dims [1], dims [2]);

    if ((chunk_dims [0] + chunk_origin [0] > dims [0]) ||
        (chunk_dims [1] + chunk_origin [1] > dims [1]) ||
        (chunk_dims [2] + chunk_origin [2] > dims [2]))
    {
      printf ("ERROR: Region overlaps with global dataset dimensions, "
          "ignoring dataset %d\n\n", i);
      continue;
    }

    data = malloc (IOnBytes (datatype, rank, chunk_dims));
    if (IOreadChunk (infile, chunk_dims, chunk_origin, data) <0)
    {
      printf ("IOreadChunk failed\n");
    }
    IOwrite (outfile, datatype, rank, chunk_dims, data);
    free (data);

    if ((j = IOreadAttributeInfo (infile, "delta", &attrType, &attrLen)) < 0)
    {
      printf ("ERROR: no delta attribute found\n\n");
    }
    else
    {
      IOreadAttribute (infile, j, &delta);
    }

    for (j = 0; j < nAttributes; j++)
    {
      if (IOreadIndexedAttributeInfo (infile, j, attrName, &attrType,
                                      &attrLen, MAXNAMESIZE) < 0)
      {
#if  0
        /* it's only a bug if this returns -1 */
        printf ("Cannot read info of attribute %d, "
            "skipping ...\n", j);
        continue;
#endif
      }
      attrData = malloc (IOnBytes (attrType, 1, &attrLen));
      if (IOreadAttribute (infile, j, attrData) < 0)
      {
        printf ("ERROR: Cannot read value of attribute %d, "
                "skipping attribute...\n", j);
        continue;
      }

      if (! strcmp (attrName, "min_ext") || ! strcmp (attrName, "origin"))
      {
        min_ext = *(double3_t *) attrData;
        min_ext.x += chunk_origin [0] * delta.x;
        min_ext.y += chunk_origin [1] * delta.y;
        min_ext.z += chunk_origin [2] * delta.z;
        *(double3_t *) attrData = min_ext;
      }
      if (! strcmp (attrName, "max_ext"))
      {
        max_ext = (double3_t *) attrData;
        max_ext->x -= (dims [0] - chunk_origin [0] - chunk_dims [0]) * delta.x;
        max_ext->y -= (dims [1] - chunk_origin [1] - chunk_dims [1]) * delta.y;
        max_ext->z -= (dims [2] - chunk_origin [2] - chunk_dims [2]) * delta.z;
      }

      IOwriteAttribute (outfile, attrName, attrType, attrLen,attrData);
      free (attrData);
    }
  }

  IOclose (infile);
  IOclose (outfile);

  return (0);
}