/***********************************************************************/
/* 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"          */
/***********************************************************************/

#include <dxconfig.h>



#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "objectClass.h"     /* for _dxfPartition */
#include "internals.h"

#define PARALLEL 1

/*
 * DXPartition - divide an a object into about n
 * spatially local partitions no smaller than size
 */

/* XXX - find a place for this */
Group _dxf_PartitionRegular(Field f, int n, int size);

Group
DXPartition(Field f, int n, int size)
{
    Group g;
    Object *o;
    int i;

    /* empty field check */
    if (DXGetObjectClass((Object)f)==CLASS_FIELD && DXEmptyField(f))
	return (Group)f /*yuck*/;

    /* try for regular field */
    g = _dxf_PartitionRegular(f, n, size);
    if (g || DXGetError()!=ERROR_NONE) {
	if (g)
	    DXCopyAttributes((Object)g, (Object)f);
	return g;
    }

    /* array to hold resulting objects */
    o = (Object *) DXAllocateZero(n * sizeof(Object));
    if (!o)
	goto error;

    /* do the partitioning */
    if (PARALLEL) {
	DXCreateTaskGroup();
	if (!_dxfPartition((Object)f, &n, size, o, 0))
	    goto error;
	if (!DXExecuteTaskGroup())
	    goto error;
    } else {
	if (!_dxfPartition((Object)f, &n, size, o, 0))
	    goto error;
    }

    /* put objects in group */
    g = (Group) DXNewCompositeField();
    for (i=0; i<n; i++)
	if (o[i])
	    DXSetMember(g, NULL, o[i]);
    DXFree((Pointer)o);
    if (g)
	DXCopyAttributes((Object)g, (Object)f);
    return g;

error:
    DXFree((Pointer)o);
    return NULL;
}


/*
 * Note - voting assumes (1) POSITIVE and NEGATIVE sum to 0 and
 * (2) all votes satisfy one and only one of ISPOS and ISNEG
 */

#define UNDEF      (-1)       /* point number not assigned yet */
#define POSITIVE   (1)        /* point is on positive side */
#define NEGATIVE   (2)        /* point is on negative side */
#define ISPOS(x)   (x>0)      /* is point on pos/neg side? */
#define ISNEG(x)   (!ISPOS(x))

#define DEP_NONE	0
#define DEP_POSITIONS   1
#define DEP_CONNECTIONS 2

#define REF_NONE	0
#define REF_POSITIONS   1
#define REF_CONNECTIONS 2

static Error
partition(Field f, int *n, int size, Object *o, int delete)
{
    struct ptinfo {		/* info about points */
	int posneg;		/* +1/-1 */
	int pos;		/* point id in pos partition */
	int neg;		/* point id in neg partition */
    } *ptinfo;
    ArrayHandle points = NULL;	/* the positions component */
    Pointer pbuf = NULL;
    ArrayHandle connections = NULL;	/* the connections component */
    Pointer cbuf = NULL;
    ArrayHandle comp = NULL;
    int *abuf = NULL;
    int npoints, nconnections;
    int npospts=0, nnegpts=0;	/* number of points in pos/neg partition */
    int nposelts=0, nnegelts=0;	/* number of connections in pos/neg partition */
    Field pos = NULL, neg = NULL;
    int i, j, k, vote, maxdim, npos, nneg, rank, dim, histo, vPerE;
    float mid, min[100], max[100], sum[100];
    Array a, pA, cA;
    char *name;
    Type type;
    Category category;
    Object attr;
    Array posA = NULL, negA = NULL;
    char *pptr, *nptr;
    int  *ipsrc = NULL, *insrc = NULL;
    byte *c_part  = NULL; /* table indicating which side the element lies in */
    int  *c_index = NULL; /* table indicating index of connection in partition */
    int  shape[32];

    /* XXX */
    histo = getenv("HISTO")? 1 : 0;

    /* get and check the positions component */
    pA = (Array) DXGetComponentValue(f, POSITIONS);
    cA = (Array) DXGetComponentValue(f, CONNECTIONS);

    if (!pA)
    {
	DXSetError(ERROR_INVALID_DATA, "field has no positions component");
	goto error;
    }

    DXGetArrayInfo(pA, &npoints, &type, &category, &rank, &dim);
    if (type!=TYPE_FLOAT || category!=CATEGORY_REAL
				   || dim<1 || dim>100 || rank != 1)
    {
	DXSetError(ERROR_INVALID_DATA, "positions have bad type");
	goto error;
    }

    if (cA)
    {
	DXGetArrayInfo(cA, &nconnections, &type, &category, &rank, &vPerE);
	if (type!=TYPE_INT || category!=CATEGORY_REAL || dim<1 || dim>100)
	{
	    DXSetError(ERROR_INVALID_DATA, "connections have bad type");
	    goto error;
	}
    }
    else
	nconnections = 0;

    if (npoints<=size || *n<=1) {
	DXEndField(f);
	*o = (Object) f;
	*n = 1;
	return OK;
    }

    points = DXCreateArrayHandle(pA);
    if (!points)
	return ERROR;

    pbuf = DXAllocateLocal(DXGetItemSize(pA));
    if (! pbuf)
	goto error;
    
    /* initialize min, max, sum */
    for (j=0; j<dim; j++) {
	min[j] = DXD_MAX_FLOAT;
	max[j] = -DXD_MAX_FLOAT;
	sum[j] = 0;
    }

    /* find bounds, sum */
    for (i=0; i<npoints; i++)
    {
	float *p = (float *)DXGetArrayEntry(points, i, pbuf);

	for (j=0; j<dim; j++)
	{
	    float m = p[j];
	    if (m < min[j]) min[j] = m;
	    if (m > max[j]) max[j] = m;
	    sum[j] += m;
	}
    }

    /* find longest dimension and its average point */
    for (maxdim=0, j=1; j<dim; j++)
	if (max[j]-min[j] > max[maxdim]-min[maxdim])
	    maxdim = j;

#define NHISTO 1000

    if (histo)
    {
	int histo[NHISTO+1], i, n;
	float scale, mn=min[maxdim], mx=max[maxdim];
	for (i=0; i<=NHISTO; i++)
	    histo[i] = 0;
	scale = NHISTO / (mx-mn);

	for (i=0; i<npoints; i++)
	{
	    float *p = (float *)DXGetArrayEntry(points, i, pbuf);
	    histo[(int)((p[maxdim]-mn)*scale)]++;
	}

	for (i=0, n=0; i<=NHISTO && n<npoints/2; i++)
	    n += histo[i];
	mid = i/scale + mn;
    }
    else
    {
	mid = sum[maxdim] / npoints;
    }

    /*
     * count pts on either side, record local point id
     */
    ptinfo = (struct ptinfo *)DXAllocateLocal(npoints*sizeof(struct ptinfo));
    if (!ptinfo) 
	goto error;

    for (i=0, nnegpts=npospts=0; i<npoints; i++)
    {
	float *p = (float *)DXGetArrayEntry(points, i, pbuf);

	if (p[maxdim] < mid)
	{
	    ptinfo[i].posneg = NEGATIVE;
	    ptinfo[i].pos = UNDEF;
	    ptinfo[i].neg = nnegpts++;
	}
	else
	{
	    ptinfo[i].posneg = POSITIVE;
	    ptinfo[i].pos = npospts++;
	    ptinfo[i].neg = UNDEF;
	}
    }

    DXFreeArrayHandle(points);
    points = NULL;
    DXFree(pbuf);
    pbuf = NULL;

    if (nconnections)
    {
	/*
	 * Go through interpolation elements assigning each to one side or
	 * the other by majority vote of which side the control points are
	 * on.  Resolve ties consistently to one partition or the other.
	 * If necessary, assign point numbers for each of the other points
	 * on the side that we just put the element on.
	 */

	c_part  = (byte *)DXAllocateLocal(nconnections*sizeof(byte));
	c_index = (int  *)DXAllocateLocal(nconnections*sizeof(int));
	if (! c_part || ! c_index)
	    goto error;

	connections = DXCreateArrayHandle(cA);
	if (! connections)
	    goto error;
    
	cbuf = (int *)DXAllocateLocal(DXGetItemSize(cA));
	if (! cbuf)
	    goto error;

	nposelts = nnegelts = 0;
	for (i = 0; i < nconnections; i++)
	{
	    int *e, *elt = (int *)DXGetArrayEntry(connections, i, cbuf);

	    for (vote=0, k=0, e=elt; k<vPerE; k++, e++)
		if (ptinfo[*e].posneg == POSITIVE)
		    vote += 1;
		else if (ptinfo[*e].posneg == NEGATIVE)
		    vote -= 1;

	    if (ISPOS(vote))
	    {
		for (k=0, e=elt; k<vPerE; k++, e++)
		{
		    struct ptinfo *p = ptinfo + *e;	
		    if (p->pos==UNDEF)		
			p->pos = npospts++;
		}

		if (c_part)
		    c_part[i] = POSITIVE;
		if (c_index)
		    c_index[i] = nposelts++;

	    }
	    else
	    {
		for (k=0, e=elt; k<vPerE; k++, e++)
		{
		    struct ptinfo *p = ptinfo + *e;
		    if (p->neg==UNDEF)		
			p->neg = nnegpts++;
		}
		if (c_part)
		    c_part[i] = NEGATIVE;
		if (c_index)
		    c_index[i] = nnegelts++;

	    }
	}

	DXFreeArrayHandle(connections);
	connections = NULL;

	DXFree(cbuf);
	cbuf = NULL;

	if (nposelts == 0 || nnegelts == 0)
	{
	    if (! DXEndField(f))
		goto error;
	    *o = (Object)f;
	    *n = 1;
	    DXFree((Pointer)c_part);
	    DXFree((Pointer)c_index);
	    return OK;
	}

    }


    /* XXX - check for no points or no faces & don't do */

    /*
     * positive, negative partitions
     * don't copy attributes: they are at top level
     * some attributes (e.g. color multipliers) must not
     * be replicated at multiple levels
     */
    pos = (Field) DXCopy((Object)f, COPY_STRUCTURE);
    while (DXGetEnumeratedAttribute((Object)pos, 0, &name))
	DXSetAttribute((Object)pos, name, NULL);
    neg = (Field) DXCopy((Object)f, COPY_STRUCTURE);
    while (DXGetEnumeratedAttribute((Object)neg, 0, &name))
	DXSetAttribute((Object)neg, name, NULL);
    
    /*
     * Now go through the components and divide them into two partitions
     */
    i = 0;
    while (NULL != (a = (Array)DXGetEnumeratedComponentValue(f, i++, &name)))
    {
	int ref = REF_NONE, dep = DEP_NONE;
	Object attr;
	int  nEntries, itemSize, nItems;

	posA = NULL;
	negA = NULL;

	if (! strcmp(name, "connections"))
	{
	    dep = DEP_CONNECTIONS;
	}
	else if (! strcmp(name, "positions"))
	{
	    dep = DEP_POSITIONS;
	}
	else
	{
	    attr = DXGetAttribute((Object)a, "dep");
	    if (attr || !strcmp(name, "connections")) 
	    {
		if (DXGetObjectClass(attr) != CLASS_STRING)
		{
		    DXSetError(ERROR_INVALID_DATA, "#10200", "dep attribute");
		    goto error;
		}

		if (! strcmp(DXGetString((String)attr), "connections"))
		    dep = DEP_CONNECTIONS;
		else if (! strcmp(DXGetString((String)attr), "positions"))
		    dep = DEP_POSITIONS;
	    }
	}
		    
	attr = DXGetAttribute((Object)a, "ref");
	if (attr)
	{
	    if (DXGetObjectClass(attr) != CLASS_STRING)
	    {
		DXSetError(ERROR_INVALID_DATA, "#10200", "ref attribute");
		goto error;
	    }

	    if (! strcmp(DXGetString((String)attr), "connections"))
		ref = REF_CONNECTIONS;
	    else if (! strcmp(DXGetString((String)attr), "positions"))
		ref = REF_POSITIONS;
	}

	if (dep == DEP_NONE && ref == REF_NONE)
	    continue;

	DXGetArrayInfo(a, &nItems, &type, &category, &rank, shape);
	itemSize = DXGetItemSize(a);
	nEntries = itemSize / DXTypeSize(type);

	/*
	 * handle constant arrays separately
	 */
	if (DXQueryConstantArray(a, NULL, NULL))
	{
	    int np, nn;

	    if (dep == DEP_POSITIONS)
	    {
		np = npospts; nn = nnegpts;
	    }
	    else
	    {
		np = nposelts; nn = nnegelts;
	    }

	    if (ref == REF_NONE)
	    {
		posA = (Array)DXNewConstantArrayV(np,
		    DXGetConstantArrayData(a), type, category, rank, shape);

		negA = (Array)DXNewConstantArrayV(nn,
		    DXGetConstantArrayData(a), type, category, rank, shape);
	    }
	    else if (ref == REF_POSITIONS)
	    {
		int *src = (int *)DXGetConstantArrayData(a);
		
		abuf = (int *)DXAllocateLocal(DXGetItemSize(a));
		if (! abuf)
		    goto error;

		for (j = 0; j < nEntries; j++)
		    if (ptinfo[src[j]].pos == UNDEF)
			abuf[j] = -1;
		    else abuf[j] = ptinfo[src[j]].pos;

		posA = (Array)DXNewConstantArrayV(np,
			    (Pointer)abuf, type, category, rank, shape);

		for (j = 0; j < nEntries; j++)
		    if (ptinfo[src[j]].neg == UNDEF)
			abuf[j] = -1;
		    else abuf[j] = ptinfo[src[j]].neg;

		negA = (Array)DXNewConstantArrayV(nn,
			    (Pointer)abuf, type, category, rank, shape);
		
		DXFree((Pointer)abuf);
		abuf = NULL;
	    }
	    else if (ref == REF_CONNECTIONS)
	    {
		int *src = (int *)DXGetConstantArrayData(a);
		
		abuf = (int *)DXAllocateLocal(DXGetItemSize(a));
		if (! abuf)
		    goto error;

		for (j = 0; j < nEntries; j++)
		    if (c_part[src[j]] == NEGATIVE)
			abuf[j] = -1;
		    else abuf[j] = c_index[src[j]];

		posA = (Array)DXNewConstantArrayV(np,
			    (Pointer)abuf, type, category, rank, shape);

		for (j = 0; j < nEntries; j++)
		    if (c_part[src[j]] == POSITIVE)
			abuf[j] = -1;
		    else abuf[j] = c_index[src[j]];

		negA = (Array)DXNewConstantArrayV(nn,
			    (Pointer)abuf, type, category, rank, shape);

		DXFree((Pointer)abuf);
		abuf = NULL;
	    }
	}
	else
	{
	    int np, nn;

	    DXGetArrayInfo(a, &nItems, &type, &category, &rank, shape);

	    comp = DXCreateArrayHandle(a);
	    if (! comp)
		goto error;

	    abuf = (int *)DXAllocateLocal(DXGetItemSize(a));
	    if (! abuf)
		goto error;

	    nn = np = 0;

	    /*
	     * Figure out how many go to either side.  If the component deps 
	     * positions or connections, we already know how many.  Otherwise,
	     * the component is referential and we have to count.  We add an
	     * element to one or both partitions depending on whether they have
	     * at least one referent into that partition.  Referents to the 
	     * other partition get a -1.
	     */
	    if (dep == DEP_POSITIONS)
	    {
		np = npospts; nn = nnegpts;
	    }
	    else if (dep == DEP_CONNECTIONS)
	    {
		np = nposelts; nn = nnegelts;
	    }
	    else if (ref == REF_POSITIONS)
	    {
		for (j = 0; j < nItems; j++)
		{
		    int p = 0, n = 0;
		    int *elt = (int *)DXGetArrayEntry(comp, j, abuf);

		    for (k = 0; k < nEntries; k++, elt++)
		    {
			if (ptinfo[*elt].pos != UNDEF) p++;
			if (ptinfo[*elt].neg != UNDEF) n++;
		    }

		    if (n) nn++;
		    if (p) np++;
		}
	    }
	    else if (ref == REF_CONNECTIONS)
	    {
		for (j = 0; j < nItems; j++)
		{
		    int p = 0, n = 0;
		    int *elt = (int *)DXGetArrayEntry(comp, j, abuf);

		    for (k = 0; k < nEntries; k++, elt++)
		    {
			if (c_part[*elt] == POSITIVE) p++;
			if (c_part[*elt] == NEGATIVE) n++;
		    }

		    if (n) nn++;
		    if (p) np++;
		}
	    }
		
	    /* 
	     * allocate the sub-arrays.
	     */
	    if (np)
	    {
		posA = (Array)DXNewArrayV(type, category, rank, shape);
		if (! posA)
		    goto error;

		if (! DXAddArrayData(posA, 0, np, NULL))
		    goto error;

		pptr = (char *)DXGetArrayData(posA);
	    }
	    else
		posA = NULL;

	    if (nn)
	    {
		negA = (Array)DXNewArrayV(type, category, rank, shape);
		if (! negA)
		    goto error;

		if (! DXAddArrayData(negA, 0, nn, NULL))
		    goto error;
		
		nptr = (char *)DXGetArrayData(negA);
	    }
	    else
		negA = NULL;

	    /*
	     * If the component contains references, we need to create buffers
	     * to re-map the references.
	     */
	    if (ref != REF_NONE)
	    {
		ipsrc = (int *)DXAllocateLocal(itemSize);
		insrc = (int *)DXAllocateLocal(itemSize);
		if (! ipsrc || ! insrc)
		    goto error;
	    }
	    else insrc = ipsrc = NULL;


	    /*
	     * Now we partition based on the dependency of the component.
	     * If the component is dep either positions or connections, we 
	     * have to put the new elements in the right place.  Otherwise,
	     * we just add 'em on to the end
	     */
	    for (j = 0; j < nItems; j++)
	    {
		char *src = (char *)DXGetArrayEntry(comp, j, abuf);
		char *psrc = (char *)src, *nsrc = (char *)src;
		int  n = 0, p = 0;

		if (ref == REF_POSITIONS)
		{
		    int *isrc = (int *)src;

		    psrc = (char *)ipsrc;
		    nsrc = (char *)insrc;

		    for (k = 0; k < nEntries; k++)
		    {
			if (isrc[k] == -1 || ptinfo[isrc[k]].pos == UNDEF)
			    ipsrc[k] = -1;
			else
			{
			    ipsrc[k] = ptinfo[isrc[k]].pos;
			    p++;
			}

			if (isrc[k] == -1 || ptinfo[isrc[k]].neg == UNDEF)
			    insrc[k] = -1;
			else
			{
			    insrc[k] = ptinfo[isrc[k]].neg;
			    n++;
			}
		    }
		}
		else if (ref == REF_CONNECTIONS)
		{
		    int *isrc = (int *)src;

		    psrc = (char *)ipsrc;
		    nsrc = (char *)insrc;

		    for (k = 0; k < nEntries; k++)
			if (isrc[k] >= 0 && c_part[isrc[k]] == POSITIVE)
			{
			    ipsrc[k] = c_index[isrc[k]];
			    insrc[k] = -1;
			    p++;
			}
			else if (isrc[k] >= 0 && c_part[isrc[k]] == NEGATIVE)
			{
			    insrc[k] = c_index[isrc[k]];
			    ipsrc[k] = -1;
			    n++;
			}
			else
			{
			    ipsrc[k] = -1;
			    insrc[k] = -1;
			}
		}
		    
		if (dep == DEP_POSITIONS)
		{
		    if (ptinfo[j].pos != UNDEF)
			memcpy(pptr+(ptinfo[j].pos*itemSize), psrc, itemSize);

		    if (ptinfo[j].neg != UNDEF)
			memcpy(nptr+(ptinfo[j].neg*itemSize), nsrc, itemSize);
		}
		else if (dep == DEP_CONNECTIONS)
		{
		    if (c_part[j] == POSITIVE)
			memcpy(pptr+(c_index[j]*itemSize), psrc, itemSize);

		    if (c_part[j] == NEGATIVE)
			memcpy(nptr+(c_index[j]*itemSize), nsrc, itemSize);
		}
		else 
		{
		    if (p)
		    {
			memcpy(pptr, psrc, itemSize);
			pptr += itemSize;
		    }

		    if (n)
		    {
			memcpy(nptr, nsrc, itemSize);
			nptr += itemSize;
		    }
		}
	    }

	    DXFreeArrayHandle(comp);

	    DXFree((Pointer)ipsrc);
	    DXFree((Pointer)insrc);
	    DXFree((Pointer)abuf);
	}

	if (! DXSetComponentValue(pos, name, (Object)posA))
	    goto error;

	if (! DXSetComponentValue(neg, name, (Object)negA))
	    goto error;
    }

    DXDeleteComponent(pos, "box");
    DXDeleteComponent(neg, "box");
    DXDeleteComponent(pos, "data statistics");
    DXDeleteComponent(neg, "data statistics");

    /* delete parent if requested, free ptinfo */
    if (delete)
	DXDelete((Object) f);

    DXFree((Pointer)ptinfo);
    DXFree((Pointer)c_part);
    DXFree((Pointer)c_index);

    /* recursively partition and delete pos */
    npos = (*n+1)/2;
    if (!_dxfPartition((Object) pos, &npos, size, o, 1))
	goto error;

    /* recursively partition and delete neg */
    nneg = *n-npos;
    if (!_dxfPartition((Object) neg, &nneg, size, o+npos, 1))
	goto error;

    /* return number of partitions generated */
    *n = npos + nneg;
    return OK;

error:

    DXFreeArrayHandle(comp);
    DXFreeArrayHandle(points);
    DXFreeArrayHandle(connections);
    DXDelete((Object)pos);
    DXDelete((Object)neg);
    DXDelete((Object)posA);
    DXDelete((Object)negA);
    DXFree(ipsrc);
    DXFree(insrc);
    DXFree(abuf);
    DXFree(pbuf);
    DXFree(cbuf);
    return ERROR;
}


/*
 * Parallelize
 */

struct arg {
    Field f;
    int n;
    int size;
    Object *o;
    int delete;
};

static Error
task(Pointer p)
{
    struct arg *arg = (struct arg *)p;
    return partition(arg->f, &arg->n, arg->size, arg->o, arg->delete);
}


Error
_dxfField_Partition(Field f, int *n, int size, Object *o, int delete)
{
    if (PARALLEL) {
	struct arg arg;
	arg.f = f;
	arg.n = *n;
	arg.size = size;
	arg.o = o;
	arg.delete = delete;
	if (!DXAddTask(task, (Pointer)&arg, sizeof(arg), 0.0))
	    return ERROR;
	return OK;
    } else {
	return partition(f, n, size, o, delete);
    }
}