/* DX module to read XmdvTool format multivariate data and create 3 dx fields
   which will form the scatterplot matrix view, the data itself, the labels,
   and the axes.  Used in scatter.net
   Last modified: 6/13/96 by MOW */

#include <stdio.h>
#include <string.h>
#include <dx/dx.h>

#define MAXDIM 20    /* arbitrary, to avoid dynamic allocation of arrays */
#define SCALE 100.0  /* this allows text to be of proper proportion to plot */


Error m_ReadMDVpc(Object *input, Object *out)
{
  FILE *in;	/* input data file */
  int i, j, k; 
  float datum, norm[MAXDIM], *p_dp, *p_lp, *p_ap;
  float llx, lly, urx, ury;
  int dim, count, *p_ac;
  char *filename, *p_ld;
  Array dp = NULL, lp = NULL, ld = NULL, ap = NULL, ac = NULL;
  Field d = NULL, l = NULL, a = NULL;
  float min[MAXDIM], max[MAXDIM], xpos[MAXDIM];
  int card[MAXDIM];
  char label[MAXDIM][80], *label_p[MAXDIM*2];
  
  if(!input[0])	{
    DXSetError(ERROR_BAD_PARAMETER, "missing filename");
    goto error;
    }
  else if(!DXExtractString(input[0], &filename))	{
    DXSetError(ERROR_BAD_PARAMETER, "filename must be a string");
    goto error;
    }
  in = fopen(filename,"r");
  if (!in) {
    fprintf(stderr, "cannot open filename %s\n", filename);
    return 0;
  }
  
/* read first line of input and insure it is valid */
  fscanf(in, "%d %d", &dim, &count);
  if (dim < 1 || dim >= MAXDIM || count < 1) {
    fprintf(stderr, "bad dimension size or count: %d, %d\n", dim, count);
    return 0;
  }
  
/* skip end of line from previous read, then scan in labels */
  fgets(label[0], 80, in); 
  for(i = 0;i < dim;i++)	{
    fgets(label[i], 80, in); 
/* strip off newline character */
    label[i][strlen(label[i])-1] = '\0';
/* store 2 links to label - one for vertical, one for horizontal */
    label_p[i*2] = label[i];
    label_p[i*2 + 1] = label[i];
    }

/* read in min and max for each dimension */
  for(i = 0;i < dim;i++)
    fscanf(in, "%f %f %d", &min[i], &max[i], &card[i]);

/* create 5 arrays */
  dp = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 2);
  lp = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 2);
  ld = DXMakeStringListV(dim*2, label_p);
  ap = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 2);
  ac = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 2);
  if(!dp || !lp || !ld || !ap || !ac)
    goto error;

/* create 3 fields */
  d = DXNewField();
  l = DXNewField();
  a = DXNewField();
  if(!d || !l || !a)
    goto error;
  
/* allocate space in arrays to hold data */
  if(!DXAddArrayData(dp, 0, dim*dim*count, NULL)) goto error;
  if(!DXAddArrayData(lp, 0, dim*2, NULL)) goto error;
  if(!DXAddArrayData(ap, 0, (dim+1)*4, NULL)) goto error;
  if(!DXAddArrayData(ac, 0, (dim+1)*2, NULL)) goto error;

/* read the data, normalize it (between 0 and 1), and write it dim**2
   times to the first output field */
  p_dp = (float *)DXGetArrayData(dp);
  if(!p_dp) goto error;

  for(i = 0;i < count;i++)	{
    for(j = 0;j < dim;j++)	{
	fscanf(in, "%f", &datum);
	norm[j] = (datum - min[j]) / (max[j] - min[j]);
	}
/* each pair of dimensions (j, k) form a plot - translate to correct place */
    for(j = 0;j < dim;j++)	{
        llx = SCALE * (float)(j) / ((float)dim);
        urx = SCALE * (float)(j+1) / ((float)dim);
        for(k = 0;k < dim;k++)	{
          lly = SCALE * (float)(k) / ((float)dim);
          ury = SCALE * (float)(k+1) / ((float)dim);
	  *(p_dp++) = llx + norm[j] * (urx - llx);
	  *(p_dp++) = lly + norm[k] * (ury - lly);
          }
        }
    }

/* finish up field description */
  if(!DXSetComponentValue(d, "positions", (Object)dp)) goto error;

  if(!DXEndField(d)) goto error;
  out[0] = (Object)d;
  
/* now create field for text labels (positions and data) */
/* note the y-value of odd position and x-value of even position is fixed */
  p_lp = (float *)DXGetArrayData(lp);
  if(!p_lp) goto error;
  for(j = 0;j < dim;j++)	{
     llx = SCALE * ((float)(j)) / ((float)dim);
     urx = SCALE * ((float)(j+1)) / ((float)dim);
     *(p_lp++) = (llx + urx) / 2.0;
     *(p_lp++) = 1.03 * SCALE; /* move it up a bit to avoid occlusion */
     *(p_lp++) = 1.03 * SCALE; /* move it over a bit to avoid occlusion */
     *(p_lp++) = (llx + urx) / 2.0;
     }


  if(!DXSetComponentValue(l, "positions", (Object)lp)) goto error;
  if(!DXSetComponentValue(l, "data", (Object)ld)) goto error;
  if(!DXEndField(l)) goto error;
  out[1] = (Object)l;
  
/* finally, create field for axes (positions and connections) */
  p_ap = (float *)DXGetArrayData(ap);
  if(!p_ap) goto error;
/* each direction has 2*(dim+1) positions and dim+1 connections */
  for(i = 0;i < dim+1;i++)	{
    llx = SCALE * ((float)(i)) / ((float)dim);
/* vertical line */
    *(p_ap++) = llx;
    *(p_ap++) = SCALE;
    *(p_ap++) = llx;
    *(p_ap++) = 0.0;
/* horizontal line */
    *(p_ap++) = SCALE;
    *(p_ap++) = llx;
    *(p_ap++) = 0.0;
    *(p_ap++) = llx;
    }

/* connections are just pair-wise */
  p_ac = (int *)DXGetArrayData(ac);
  if(!p_ac) goto error;
  for(i = 0;i < (dim+1)*4;i= i+2)	{
    *(p_ac++) = i;
    *(p_ac++) = i+1;
    }

  if(!DXSetStringAttribute((Object)ac, "element type", "lines")) goto error;
  if(!DXSetComponentValue(a, "positions", (Object)ap)) goto error;
  if(!DXSetComponentValue(a, "connections", (Object)ac)) goto error;
  if(!DXEndField(a)) goto error;
  out[2] = (Object)a;

  return OK;
/* clean up everything on error */
error:
  DXDelete((Object)d);
  DXDelete((Object)l);
  DXDelete((Object)a);
  DXDelete((Object)dp);
  DXDelete((Object)lp);
  DXDelete((Object)ld);
  DXDelete((Object)ap);
  DXDelete((Object)ac);
}