/*
 *
 * Convert MPP format movie and history dumps
 * to something that micomproc can swallow.
 *
 * Compliation instructions: "cc -O -o cvtmpp cvtmpp.c"
 *
 * Aaron Sawdey, 1996
 *
 */

#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>

#if defined(_CRAY) && !defined(_CRAYMPP)
#define VECTOR
#endif

#define MALLOC(size) chk_malloc(size,__LINE__,__FILE__)

extern char *optarg;
extern int optind;
int debug=0;

char *progname;
char *outputdir = ".";
int altpakk=0;		/* set to 1 for "linda's" pakk format */
int cvtpakk=0;

int idm,jdm;

/*
 * All the info we might need about
 * a processors piece of data.
 */
struct piece_info {
  char cname[10];	/* field name */
  int pe;		/* processor # */
  long nstep;		/* step number */
  double time,coord;	/* time and sigma coord */
  int is,js,ie,je;	/* location */
  int len;		/* length of packed data */
  int fno;		/* file number */
  off_t offset;		/* offset in file */
};

struct sect_info {
  int is,js,ie,je;
};

struct piece_info *readpiece(FILE *, char *, char *);
  
/*
 * fprintf(coutf," ::%8s%4d%8d%9.2f%8.5f%4d%4d%4d%4d%6d\n",
 *	   cname,who,nstep,time,coord,is,ie,js,je,len);
 */

void
usage()
{
  fprintf(stderr,
	  "Usage: %s [-d] [-s idm,jdm] files\n"
	  "  -a             expect 'alternate' pakk data format.\n"
	  "  -c		    convert between the two pakk formats.\n"
	  "                 (one is used for input, the other for output)\n"
	  "  -d             turn on debug code.\n"
	  "  -n npes        number of processors used.\n"
	  "  -o dir         directory to put output files in.\n"
	  "  -s idm,jdm     set size of dataset.\n"
	  "  -S si,sj       split output into si by sj pieces.\n"
	  "  -P is,js,ie,je output one piece containing the region\n"
	  "                 (is,js)-(ie,je).\n",
	  progname);
  

  exit(1);
}

void *
chk_malloc(long size, int line, char *file)
{
  void *rptr;

  rptr = (void*)malloc(size);

  if(rptr == NULL) 
    {
      fprintf(stderr,
	      "Out of memory allocating %d bytes at line %d in file %s.\n",
	      size,line,file);
      exit(1);
    }	
  return rptr;
}

/*
 * Produce an array of sect_info structures describing
 * the desired geometry of pieces of ocean.
 */
struct sect_info *
makesects(int pieces_i, int pieces_j)
{
  struct sect_info *si;
  int i,j,k;

  si = MALLOC(pieces_i*pieces_j*sizeof(struct sect_info));

  k = 0;
  for(j=0;j<pieces_j;j++)
    for(i=0;i<pieces_i;i++)
      {
	si[k].is = (i*idm)/pieces_i;
	si[k].ie = ((i+1)*idm)/pieces_i - 1;
	si[k].js = (j*jdm)/pieces_j;
	si[k].je = ((j+1)*jdm)/pieces_j - 1;
	k++;
      }
  return si;
}

/*
 * Some simple routines for parsing things written
 * by fortran that may have fields running together (which
 * sscanf doesn't handle well at all).
 */

static char *f_ptr;
void f_new(char *str) {f_ptr = str; return;}
void f_skip(int n) {f_ptr += n;}
void f_char(char *buf, int n) {int i; for(i=0;i<n;i++) buf[i] = f_ptr[i];
			       buf[n] = '\0'; f_ptr += n;}
int f_int(int n) {char buf[100]; f_char(buf,n); return(atoi(buf));}
double f_real(int n) {char buf[100]; f_char(buf,n); return(atof(buf));}

struct piece_info *dir;	/* directory of pieces. */
int nfiles;
char **fnames;
FILE **fps;
int npes = 0;
int fields;
int maxpi,maxpj; /* max size of piece */
off_t *fsect;
off_t *flen;
float *array;
char *piecebuf;

#ifdef VECTOR
typedef long iarray_t;
#else
typedef int iarray_t;
#endif
iarray_t *iarray;

int
main(int argc, char **argv)
{
  int c;
  int i,j;
  char *cp;
  int perfile;
  int done;
  struct sect_info *sects = NULL;
  int nsects;
  int pi,pj;

  progname = argv[0];
  npes = 0;

  pi = 5; pj = 5;

  if(argc == 1) usage();

  while((c = getopt(argc,argv,"acn:o:s:dS:P:")) != -1)
    {
      switch(c) {
      case 'a':	altpakk=1; break;
      case 'c':	cvtpakk=1; break;
      case 'n':
	npes = atoi(optarg);
	break;

      case 'o':
	outputdir = optarg;
	break;

      case 'S':
	{ int si,sj;
	  i = sscanf(optarg,"%d,%d",&si,&sj);
	  if(i != 2) {
	    fprintf(stderr,"The -S option takes two comma-separated "
		    "integer arguments.\n");
	    usage();
	  }
	  pi = si; pj = sj;
	  break;
	}

      case 'P':
	{ int is,js,ie,je;
	  i = sscanf(optarg,"%d,%d,%d,%d",&is,&js,&ie,&je);
	  if(i != 4) {
	    fprintf(stderr,"The -P option takes four comma-separated "
		    "integer arguments.\n");
	    usage();
	  }
	  nsects = 1;
	  sects = malloc(sizeof(struct sect_info));
	  sects->is = is-1;
	  sects->js = js-1;
	  sects->ie = ie-1;
	  sects->je = je-1;
	  break;
	}

      case 's': 
	i = sscanf(optarg,"%d,%d",&idm,&jdm);
	if(i != 2) {
	  fprintf(stderr,"The -s option takes two comma-separated "
		  "integer arguments.\n");
	  usage();
	}
	if(debug) printf("idm,jdm set to %d,%d\n",idm,jdm);
	break;

      case 'd': 
	debug++; break;

      default:
	usage();
      }
    }

  nfiles = argc-optind;

  if(nfiles == 0) usage();

  if(nfiles == 1 && npes == 0) {
    fprintf(stderr,"If there is only one file, you have to tell me\n"
	    "how many processors were used. Use the -n option.\n");
    usage();
  }

  if(debug) {
    printf("nfiles = %d\n",nfiles);
  }

  /* Allocate array of filenames and file pointers */
  fnames = (char **)MALLOC(nfiles * sizeof(char*));
  fps = (FILE **)MALLOC(nfiles * sizeof(FILE*));

  /* open all the files. */
  for(i = 0 ; i < nfiles ; i++) {
    fnames[i] = strdup(argv[optind+i]);
    fps[i] = fopen(fnames[i],"r");
    if(debug) printf("Opened '%s'.\n",fnames[i]);
    if(fps[i] == NULL) {
      fprintf(stderr,"Can't open input file '%s':\n",fnames[i]);
      perror(fnames[i]);
      exit(1);
    }
  }

  if(npes == 0) {
    cp = strrchr(fnames[0],'/');
    if(cp == NULL) cp = fnames[0];
    i = atoi(cp+1);
    cp = strrchr(fnames[1],'/');
    if(cp == NULL) cp = fnames[1];
    j = atoi(cp+1);
    i -= j;
    if(i < 0) i = -i;

    npes = i*nfiles;
    if(debug) printf("Figured out that npes = %d\n",npes);
  }

  perfile = npes/nfiles;

  dir = (struct piece_info *)MALLOC(sizeof(struct piece_info)*npes);
  fsect = (off_t*)MALLOC(sizeof(off_t)*nfiles);
  flen = (off_t*)MALLOC(sizeof(off_t)*nfiles);

  for(j=0;j<nfiles;j++)
    for(i=0;i<perfile;i++)
      {
	dir[i+j*perfile] = *readpiece(fps[j],fnames[j],NULL);
	dir[i+j*perfile].fno = j;
      }

  /* Record the length of a section in each file.
   * This should be the same for each section within
   * a particular file.
   */
  for(j=0;j<nfiles;j++)
    {
      fsect[j] = ftell(fps[j]);
      fseek(fps[j],0,2);
      flen[j] = ftell(fps[j]);
      fseek(fps[j],0,0);
    }

  idm = 0; jdm = 0;
  maxpi = 0; maxpj = 0;
  for(i=0;i<npes;i++) {
    int pi = dir[i].ie - dir[i].is + 1;
    int pj = dir[i].je - dir[i].js + 1;

    if(dir[i].ie > idm) idm = dir[i].ie;
    if(dir[i].je > jdm) jdm = dir[i].je;
    if(pi > maxpi) maxpi = pi;
    if(pj > maxpj) maxpj = pj;
  }

  if(sects == NULL) {
    nsects = pi*pj;
    sects = makesects(pi,pj);
  }

  array = (float *)MALLOC(idm*jdm*sizeof(float));
  iarray = (iarray_t *)MALLOC(2*idm*sizeof(iarray_t));
  piecebuf = (char *)MALLOC(2*idm*jdm+50+(idm*jdm)/64);

  fields = flen[0] / fsect[0];
  for(i = 0 ; i < nfiles ; i++) {
    if(fields*fsect[i] != flen[i]) {
      fprintf(stderr,"File %s is the wrong length (was %d, should be %d).\n",
	      fnames[i],flen[i],fields*fsect[i]);
      fprintf(stderr,"fields=%d flen=%d fsect=%d\n",fields,flen[i],fsect[i]);
      exit(1);
    }
  }

  if(debug) {
    printf("Allocated %d bytes for %d,%d array.\n",
	   idm*jdm*sizeof(float),idm,jdm);
    printf("Max piece is %d,%d.\n",maxpi,maxpj);
    printf("The files contain %d fields.\n",fields);
  }

  for(i=0;i<nsects;i++) {
    char on[100];
    sprintf(on,"%s/%04d,%04d",outputdir,sects[i].is+1,sects[i].js+1);
    write_sect(on,&(sects[i]));
  }

  exit(0);
}

/*
 * Return nonzero if the domains represented by p and s
 * overlap, zero otherwise.
 */
int
intersects(struct piece_info *p, struct sect_info *s)
{
  return (((p->is >= (s->is+1) && p->is <= (s->ie+1)) || 
	   (p->ie >= (s->is+1) && p->ie <= (s->ie+1))) &&
	  ((p->js >= (s->js+1) && p->js <= (s->je+1)) || 
	   (p->je >= (s->js+1) && p->je <= (s->je+1))));
}


/* 
 * This reads in the neccessary processor sections to write out 
 * the given section of the data to the specified filename.
 */
int
write_sect(char *fname, struct sect_info *sect)
{
  FILE *out;
  int f,p,n,i,j,jidm;
  char *pbuf;
  off_t o;
  struct piece_info pi;
  int lenlen=6;

  pbuf = (char *)MALLOC(maxpi*maxpj*2+1000);

  out = fopen(fname,"w");
  if(out == NULL) {
    fprintf(stderr,"Can't open output file %s.\n",fname);
    perror(fname);
    exit(1);
  }

  printf("Collecting section for %d,%d - %d,%d\n",sect->is+1,sect->js+1,
	 sect->ie+1,sect->je+1);

  for(f = 0 ; f < fields ; f++) {

    /*
     * Fill up the array with zeros.
     * Perhaps it might be more desireable to use
     * another value for some fields.
     */
    for(j=sect->js ; j <= sect->je ; j++) {
      jidm = j*idm;
      for(i=sect->is ; i <= sect->ie ; i++)
	array[i+jidm] = 0.0;
    }
      
    for(p = 0 ; p < npes ; p++) {
      if(intersects(&(dir[p]),sect)) {
	if(debug>2) {
	  printf("(%d,%d)-(%d,%d) overlaps\n",
		 dir[p].is,dir[p].js,dir[p].ie,dir[p].je);
	}
	n = dir[p].fno;
	o = f*fsect[n] + dir[p].offset;
	fseek(fps[n],o,0);
	pi = *readpiece(fps[n],fnames[n],pbuf);
	unpak(array, pbuf, &(dir[p]));
      }
    }
    i = sect->ie-sect->is+1;
    j = sect->je-sect->js+1;
    if(i*j+14 >= 1000000) lenlen = 8;

    fprintf(out," ::%-8s%8d%8.1f%8.5f%4d%4d%*d\n",
	    pi.cname,pi.nstep,pi.time,pi.coord,i,j,lenlen,i*j+14);
	    
    pak_out(out,array,sect);
  }

  free(pbuf);

  fclose(out);
  return 0;
}

/*
 * fprintf/fwrite stmts used in csupport.c:
 *
 * fprintf(coutf," ::%8s%4d%8d%9.2f%8.5f%4d%4d%4d%4d%6d\n",
 *	   cname,who,nstep,time,coord,is,ie,js,je,len);
 * fwrite(buf,2,len,coutf);
 *
 */
struct piece_info *
readpiece(FILE *fp, char *fname, char *pdata)
{
  static struct piece_info rval;
  char lbuf[1024];
  char tc[10];
  int i;
  int datasz;

  rval.offset = ftell(fp);

  if((fgets(lbuf,sizeof(lbuf),fp)) == NULL) {
    return NULL;
  }

  if(debug > 2) {
    printf("Reading from file %s, fp=0x%08x\n",fname,fp);
  }

  f_new(lbuf);
  f_char(tc,3);
  if(strcmp(tc," ::")) {
    fprintf(stderr,"Error reading file %s, expected ' ::' and got:\n%s",
	    fname,lbuf);
    exit(1);
  }

  f_char(rval.cname,8);
  rval.pe = f_int(4);
  rval.nstep = f_int(8);
  rval.time = f_real(9);
  rval.coord = f_real(8);
  rval.is = f_int(4);
  rval.ie = f_int(4);
  rval.js = f_int(4);
  rval.je = f_int(4);
  rval.len = f_int(6);

  /* number of bytes to read or skip */
  datasz = 1 + 2*rval.len;

  if(debug > 2) {
    printf("datasz = %d\n",datasz);
  }

  if(pdata == NULL) {
    fseek(fp,datasz,1);
  } else {
    i = fread(pdata,1,datasz,fp);
    if(i != datasz) {
      perror(fname);
      fprintf(stderr,"Partial read on %s.\n",fname);
      exit(1);
    }
  }

  return &rval;
}

int
unpak(float *data, char *buf, struct piece_info *pi)
{
  int i,j,l,m;
  int i1,i2;
  char baseb[16],scalb[16];
  double base,scal,fval;
  int trans[129];

  memcpy(baseb,buf,14);
  baseb[14] = '\0';

  memcpy(scalb,buf+14,14);
  scalb[14] = '\0';

#ifndef VECTOR
  for(i = 0 ; i < 129 ; i++) {
    i1 = i;
    if(!altpakk) {
      if (i1 > 96) i1=i1-6;
      if (i1 > 64) i1=i1-7;
      i1 = i1 - 14;
    }
    i1=i1-32;
    trans[i] = i1;
  }
#endif

  for(i=0;i<14;i++) {
    if(baseb[i] == 'E') baseb[i] = 'e';
    if(scalb[i] == 'E') scalb[i] = 'e';
  }

  base = atof(baseb);
  scal = atof(scalb);

  l=28;
    for(i = pi->is-1 ; i < pi->ie ; i++)
    {
#ifdef VECTOR
      m = pi->je - pi->js + 1;
      _unpack(buf+l,iarray,m*2,-1);
      l += m*2;

      for(j = pi->js-1,m=0 ; j < pi->je ; j++,m+=2)
	{
	  i1 = iarray[m];
	  i2 = iarray[m+1];
	  if(!altpakk) {
	    if (i1 > 96) i1=i1-6;
	    if (i1 > 64) i1=i1-7;
	    if (i2 > 96) i2=i2-6;
	    if (i2 > 64) i2=i2-7;
	    i1 = i1 - 14;
	    i2 = i2 - 14;
	  }
	  i1=i1-32;
	  i2=i2-32;
	  data[i+j*idm]=scal*(float)((i1<<6)|i2)+base;
	}
#else
      for(j = pi->js-1 ; j < pi->je ; j++)
	{
	  i1 = trans[buf[l++]];
	  i2 = trans[buf[l++]];
	  
	  data[i+j*idm]=scal*(float)((i1<<6)|i2)+base;
	}
#endif
    }
  return 0;
}

#define LINELEN 128

pak_out(FILE *fp, float *data, struct sect_info *si)
{
  int i,j,l,m;
  int c1,c2,c3;
  int jidm;
  int numb,i1,i2;
  int outfmt;
  int trans[129];

  float base,maxv,scal,t;
  double iscal;

  base = data[si->is + (si->js)*idm];
  maxv = base;

  outfmt = altpakk;
  if(cvtpakk) {
    if(outfmt) outfmt = 0;
    else outfmt = 1;
  }

#ifndef VECTOR
  for(i = 0 ; i < 129 ; i++) {	
    i1 = i;
    if(!outfmt) {	
      if (i1 > 37) i1=i1+6;
      if (i1 > 11) i1=i1+7;
      i1 = i1 + 14;
    }
    i1=i1+32;
    trans[i] = i1;
  }
#endif

  for(j = si->js ; j <= si->je ; j++)
    {
      jidm = j*idm;
      for(i = si->is ; i <= si->ie ; i++)
	{	
	  t = data[i + jidm];
	  if(t > maxv) maxv = t;
	  if(t < base) base = t;
	}
    }

  scal = (maxv-base)/4095.0;

  if(scal == 0.0) scal = 1.0;

  iscal = 1.0/scal;

  sprintf(piecebuf,"%14.7e%14.7e",base,scal);
  l = 28;

  for(i = si->is ; i <= si->ie ; i++)
    {
#ifdef VECTOR
      if(!outfmt) {
	c1 = 6;
	c2 = 7;
	c3 = 14;
      } else {
	c1 = 0; 
	c2 = 0;
	c3 = 0;
      }
      for(j = si->js,m=0 ; j <= si->je ; j++,m+=2)
	{
	  numb=(data[i+j*idm]-base)/scal+.5;
	  i1 = (numb>>6);
	  i2 = (numb&63);
	  if (i1 > 37) i1=i1+c1;
	  if (i1 > 11) i1=i1+c2;
	  if (i2 > 37) i2=i2+c1;
	  if (i2 > 11) i2=i2+c2;
	  i1=i1+32+c3;
	  i2=i2+32+c3;
	  iarray[m] = i1;
	  iarray[m+1] = i2;
	}

      m = si->je - si->js + 1;
      _pack(iarray,piecebuf+l,2*m,-1);
      l+=2*m;
#else
      for(j = si->js ; j <= si->je ; j++)
	{
	  numb=(data[i+j*idm]-base)/scal+.5;
	  i1 = trans[(numb>>6)];
	  i2 = trans[(numb&63)];
	  piecebuf[l++] = i1;
	  piecebuf[l++] = i2;
	}
#endif
    }

  for(i = 0 ; i < l ; i+= 128) {
    if(l-i > 128) j = 128;
    else j = l-i;
		   
    fputc(' ',fp);
    fwrite(piecebuf+i,j,1,fp);
    fputc('\n',fp);
  }

  return 0;
}