/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/
/*
 * $Header: /home/gda/dxcvs/dx/src/exec/dxmods/transpose.c,v 1.3 1999/05/10 15:45:32 gda Exp $
 */

#include <dxconfig.h>



#include <dx/dx.h>
#include <ctype.h>

static Object Object_Transpose(Object o, int num, int *dimlist);
static Field Field_Transpose(Field f, int num, int *dimlist);
static Array Item_Transpose(Array a, int num, int *dimlist, char *comp);
static Array Connect_Transpose(Array a, int num, char *element);


#define STR(x)    (Object)DXNewString(x)

#define MAXRANK   16
#define MAXSHAPE 128



#define PARALLEL 1

 
#if PARALLEL
static Field Transpose_Wrapper(Field, char *, int);
struct argblk {
    int ndim;
    int dim[MAXRANK];
};
#endif

#define is_empty(o)  ((DXGetObjectClass(o)==CLASS_FIELD) && DXEmptyField((Field)o))


/* 0 = input
 * 1 = new dimension order
 *
 * 0 = output
 */
int
m_Transpose(Object *in, Object *out)
{
    int i, ndim;
    char cdim, *cp;
    int dim[MAXRANK];

    if(!in[0]) {
	DXSetError(ERROR_BAD_PARAMETER, "#10000", "input");
	return ERROR;
    }

    if(!in[1] || is_empty(in[0])) {
	out[0] = in[0];
	return OK;
    }

    /* input can either be a list of integers, or a list of strings.
     */
    if (DXQueryParameter(in[1], TYPE_STRING, 1, &ndim)) {
	
	/* list of strings - since there isn't a get-number-of-strings call,
	 *  check first to be sure there aren't too many, and then loop 
	 *  until there aren't any more.
	 */
	if(DXExtractNthString(in[1], MAXRANK, &cp)) {
	    DXSetError(ERROR_BAD_PARAMETER, "bad number of dimensions in list");
	    return ERROR;
	}
	for(i=0; i<MAXRANK; i++) {
	    if(!DXExtractNthString(in[1], i, &cp))
		break;
	    else {
		cdim = isupper(cp[0]) ? tolower(cp[0]) : cp[0];
		if(cdim < 'x' || cdim > 'z') {
		    DXSetError(ERROR_BAD_PARAMETER, "bad dimension letter");
		    return ERROR;
		}
		dim[i] = (int)(cdim - 'x');
	    }
	}
	ndim = i;
	if(ndim <= 1  || ndim > MAXRANK) {
	    DXSetError(ERROR_BAD_PARAMETER, "bad number of dimensions in list");
	    return ERROR;
	}
	
    } else if(DXQueryParameter(in[1], TYPE_INT, 1, &ndim)) {
	
	/* list of ints */
	if(ndim <= 1  || ndim > MAXRANK) {
	    DXSetError(ERROR_BAD_PARAMETER, "bad number of dimensions in list");
	    return ERROR;
	}
	if(!DXExtractParameter(in[1], TYPE_INT, 1, ndim, (Pointer)dim)) {
	    DXSetError(ERROR_BAD_PARAMETER, "bad dimension list");
	    return ERROR;
	}
	
    } else {
	/* neither integers or strings */
	DXSetError(ERROR_BAD_PARAMETER, "bad dimension list");
	return ERROR;
    }



#if PARALLEL
    { struct argblk argblk;
      int i;
      
      argblk.ndim = ndim;
      for (i=0; i<ndim; i++)
	  argblk.dim[i] = dim[i];

      out[0] = DXProcessParts(in[0], Transpose_Wrapper, (char *)&argblk, 
			    sizeof(argblk), 1, 1);
    
    }
#else

    if (!(out[0] = DXCopy(in[0], COPY_STRUCTURE))) {
	DXSetError(ERROR_BAD_PARAMETER, "bad input object");
	return ERROR;
    }
	    
    out[0] = Object_Transpose(out[0], ndim, dim);
#endif
    
    return(out[0] ? OK : ERROR);


}

#if PARALLEL
/*
 */ 
static Field Transpose_Wrapper(Field f, char *argblk, int size)
{
    return (Field)Object_Transpose((Object)f, 
				   ((struct argblk *)argblk)->ndim,
				   ((struct argblk *)argblk)->dim);
}
#endif

/*
 */
static Object Object_Transpose(Object o, int num, int *dims)
{
    Object subo;
    char *name;
    int i;

    switch(DXGetObjectClass(o)) {
      case CLASS_GROUP:
	/* for each member */
	for(i=0; subo = DXGetEnumeratedMember((Group)o, i, &name); i++) {
	    if(!Object_Transpose(subo, num, dims))
		return NULL;
	}
	break;

      case CLASS_FIELD:
#if ! PARALLEL
        if(DXEmptyField((Field)o))
            break;

#else
        if(DXEmptyField((Field)o))
            return NULL;

	o = DXCopy(o, COPY_STRUCTURE);
	if(!o)
	    return NULL;
#endif
	if (!Field_Transpose((Field)o, num, dims)) {
	    DXDelete(o);
	    o = NULL;
	}
	break;
	
      default:
        /* ignore things we don't understand */
        break;
    }
    
    return o;
}




/*
 */
static Field Field_Transpose(Field f, int num, int *dims)
{
    Object newo, subo;
    char *name, *attr;
    int i, j, n, nitems, rank;
    int inverted = 0;
    int ndims[MAXRANK];
    int pcounts[MAXRANK], ccounts[MAXRANK];
    float origins[MAXRANK], norigins[MAXRANK];
    float deltas[MAXRANK*MAXRANK], ndeltas[MAXRANK*MAXRANK];


    /* transpose positions.
     */
    subo = DXGetComponentValue(f, "positions");
    if(!subo) {
        DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
	return ERROR;
    }
    
    if(!DXTypeCheckV((Array)subo, TYPE_FLOAT, CATEGORY_REAL, 1, NULL)) {
        DXSetError(ERROR_INVALID_DATA, "positions must be type float, rank 1"); 
        return NULL;
    }
    
    /* make sure all the dims are valid before proceeding */
    DXGetArrayInfo((Array)subo, &nitems, NULL, NULL, &rank, pcounts);
    if(rank != 1) {
        if(rank == 0)
            pcounts[0] = 1;    /* rank 0 == rank 1, shape 1 */
        else
            DXErrorReturn(ERROR_INVALID_DATA, "bad position component");
    }
    
    if(num != pcounts[0]) {
        DXSetError(ERROR_BAD_PARAMETER, "#11365", num, pcounts[0]);
        return NULL;
    }
    
    /* initialize a use-counts list.  in the second loop, if we find the 
     *  same dimension more than once in the list, it's an error.
     */
    for(i=0; i<pcounts[0]; i++)
        ndims[i] = 0;
    
    for(i=0; i<pcounts[0]; i++) {
        /* out of range?
	 */
        if(dims[i] < 0  || dims[i] >= num) {
            DXSetError(ERROR_BAD_PARAMETER, "#11361", dims[i], "input", num-1);
            return NULL;
        }
        /* duplicates?
	 */
        if(++ndims[dims[i]] > 1) {
            DXSetError(ERROR_BAD_PARAMETER, 
                     "dimension %d appears more than once in list",
                     dims[i]);
            return NULL;
        }
    }
    
    
    /* dim list should be ok.  now rearrange positions 
     */
    if(DXQueryGridPositions((Array)subo, &n, pcounts, origins, deltas)) {
        /* regular grid 
	 */
        for(i=0; i<n; i++)
            norigins[i] = origins[dims[i]];
	
        for(i=0; i<n; i++)
            for(j=0; j<n; j++)
                ndeltas[i*n + j] = deltas[i*n + dims[j]];

	
        newo = (Object)DXMakeGridPositionsV(n, pcounts, norigins, ndeltas);
        

    } else {
        /* irregular positions, have to reorder them inside each point.
	 */
        newo = (Object)Item_Transpose((Array)subo, num, dims, "positions");
        
        inverted = 1;
        for(i=0; i<num; i++)
            if(dims[i] != i)
                inverted = !inverted;
        
    }

    if (!newo)
	return ERROR;

    if (!DXSetComponentValue(f, "positions", (Object)newo)
	|| !DXChangedComponentValues(f, "positions"))
	return ERROR;
    
    
    /* connections, if they exist and need to be inverted.
     */
    subo = DXGetComponentValue(f, "connections");
    if (subo && inverted) {
	
	if (!DXTypeCheckV((Array)subo, TYPE_INT, CATEGORY_REAL, 1, NULL)) {
	    DXSetError(ERROR_INVALID_DATA, "bad connections component"); 
	    return NULL;
	}
	
	/* if the connections are irregular, and we've inverted the axes,
	 *  then we have to reorder the verticies of the connections.
	 */
	if (!DXQueryGridConnections((Array)subo, &n, ccounts)) {
	    
	    DXGetArrayInfo((Array)subo, &nitems, NULL, NULL, &rank, ccounts);
	    if (rank != 0  &&  rank != 1) {
		DXSetError(ERROR_INVALID_DATA, "connections component not rank 1");
		return NULL;
	    }
	    
	    attr = DXGetString((String)DXGetAttribute((Object)subo, "element type"));
	    if (!attr) {
		DXSetError(ERROR_INVALID_DATA, "connections have no element type");
		return NULL;
	    }
	    
	    newo = (Object)Connect_Transpose((Array)subo, ccounts[0], attr);
	    if (!newo)
		return NULL;
	    
	
	    if (!DXSetComponentValue(f, "connections", (Object)newo) ||
		!DXChangedComponentValues(f, "connections"))
		return ERROR;

	}
    }
	

    /* for each of the rest of the components, see if they need 
     *  to be transposed.
     */
    for(i=0; subo = DXGetEnumeratedComponentValue(f, i, &name); i++) {

        /* skip positions & connections.
	 */
        if (!strcmp("positions", name) || !strcmp("connections", name))
            continue;
        
	/* for now, special case normals.  as soon as we can agree,
	 *  skip any component which doesn't have the 'XXX' attribute set.
	 */
        attr = DXGetString((String)DXGetAttribute(subo, "geometric"));
	
	if (!attr && strcmp(name, "normals"))
	    continue;

	newo = (Object)Item_Transpose((Array)subo, num, dims, name);

	if (!newo ||
	    !DXSetComponentValue(f, name, (Object)newo) ||
	    !DXCopyAttributes(subo, newo))
	    return ERROR;
	
    }
    
    DXEndField(f);
    return f;

}


/*
 * reorder the connections in an array, if necessary
 */
static Array Connect_Transpose(Array a, int num, char *element)
{
    Array na = a;
    int i, ndims[MAXRANK];

    if (!strcmp(element, "triangles") || !strcmp(element, "tetrahedra")) {
	
	/* construct a new dim array with 0<->1, and call item_tr() */
	ndims[0] = 1;
	ndims[1] = 0;
	for(i=2; i<num; i++)
	    ndims[i] = i;
	
	na = Item_Transpose(a, num, ndims, "connections");
	
    } 
    else if (!strcmp(element, "cubes"))
	;
    else if (!strcmp(element, "quads"))
	;
    else if (!strcmp(element, "lines"))
	;
    else if (!strcmp(element, "points"))
	;
    else {
	DXWarning("unrecognized element type; connections not inverted");
    }

    return na;
}


/* 
 * for each vector item in an array, reorder the items inside the vector
 */
static Array Item_Transpose(Array a, int num, int *dims, char *comp)
{
    Type t;
    Category c;
    int rank, shape[MAXRANK];
    float origins[MAXRANK], norigins[MAXRANK];
    float deltas[MAXRANK], ndeltas[MAXRANK];
    Array na;
    Array terms[MAXRANK], nterms[MAXRANK];
    int nitems, typesize;
    Pointer op, np;
    float *ofp, *nfp;
    int i, j;
    

    /* regardless of the array type, these checks apply */
    if(!DXGetArrayInfo(a, &nitems, &t, &c, &rank, shape))
	return NULL;
    
    if(rank == 0 || (rank == 1 && shape[0] == 1))
	return a;        /* transpose of single wide vector is itself */
    
    if(rank != 1) {
	DXSetError(ERROR_INVALID_DATA, "can only transpose vector data");
	return NULL;     /* only transposes rank 1 arrays */
    }
    
    if(num != shape[0]) {
	DXSetError(ERROR_INVALID_DATA, 
		   "number of dimensions doesn't match shape of '%s'", comp);
	return NULL;
    }
	

    switch(DXGetArrayClass(a)) {
	
      default:
      case CLASS_ARRAY:
	/* normal expanded array */
	
	na = DXNewArrayV(t, c, rank, shape);
        if(!na)
	    return NULL;
	
	np = DXGetArrayData(DXAddArrayData(na, 0, nitems, NULL));
	op = DXGetArrayData(a);

	nitems *= DXCategorySize(c);

	/* for speed, special case floats */
	if(t == TYPE_FLOAT && c == CATEGORY_REAL) {
	    ofp = (float *)op;
	    nfp = (float *)np;

	    for(i=0; i<nitems; i++) {
		for(j=0; j<shape[0]; j++)
		    nfp[j] = ofp[dims[j]];

		nfp += shape[0];
		ofp += shape[0];
	    }
	} else {
	    typesize = DXTypeSize(t);
	    for(i=0; i<nitems; i++) {

		for(j=0; j<shape[0]; j++)
		    memcpy((char *)np + dims[j]*typesize, 
			   (char *)op + j*typesize,
			   typesize);

		np = (Pointer)((char *)np + typesize*shape[0]);
		op = (Pointer)((char *)op + typesize*shape[0]);
	    }
	}

	break;

      case CLASS_CONSTANTARRAY:

	np = DXAllocate(DXGetItemSize(a));
	op = DXGetArrayData(a);
	if (!np || !op)
	    return NULL;

	typesize = DXTypeSize(t) * DXCategorySize(c);
	for(j=0; j<shape[0]; j++)
	    memcpy((char *)np + dims[j]*typesize, 
		   (char *)op + j*typesize,
		   typesize);
	    

	na = (Array) DXNewConstantArrayV(nitems, np, t, c, rank, shape);
	DXFree(np);
	break;

      case CLASS_REGULARARRAY:

	if(!DXGetRegularArrayInfo((RegularArray)a, &i, 
				(Pointer)origins, (Pointer)deltas))
	    return NULL;

	for(j=0; j<shape[0]; j++) {
	    norigins[j] = origins[dims[j]];
	    ndeltas[j] = deltas[dims[j]];
	}
	na = (Array)DXNewRegularArray(t, shape[0], i, 
				    (Pointer)norigins, (Pointer)ndeltas);
	break;

      case CLASS_PATHARRAY:
	/* doesn't change */

	na = a;
	break;

      case CLASS_PRODUCTARRAY:
	/* recurse until you get to a more fundamental type of array */

	DXGetProductArrayInfo((ProductArray)a, &j, terms);
	for(i=0; i<j; i++)
	    nterms[i] = Item_Transpose(terms[i], num, dims, comp);

	na = (Array)DXNewProductArrayV(j, nterms);
	break;
	
      case CLASS_MESHARRAY:
	/* recurse until you get to a more fundamental type of array */

	DXGetMeshArrayInfo((MeshArray)a, &j, terms);
	for(i=0; i<j; i++)
	    nterms[i] = Item_Transpose(terms[i], num, dims, comp);

	na = (Array)DXNewMeshArrayV(j, nterms);
	break;
	
    }
    
    return na;
    
}