/***********************************************************************/ /* 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 #include #include #include /* internal statistics structure. this is what the "xxx statistics" * component actually contains, but users have to go through the * DXStatistics() call so they never reference it directly, and it's * derived (has a 'der' attribute) so it doesn't get exported, so it is * reasonably safe to change it if necessary. */ struct stats { double min, max, avg, sigma; /* the values needed to combine stats */ double sum, sumsq, counts; /* from different fields */ double isvalid; /* set to 0 if all points invalid */ }; /* to handle composite fields, where you need to know the stats of this * one piece by itself, but you also need to know the stats if this is * considered as just one part of the whole (it will be different because * duplicate points will be ignored in one and not in the other). * * the has_cf flag isn't used right now because we can tell the difference * between a component which just has stats and one which has both types * of stats by the number of items in the array. */ struct statcomp { struct stats f; /* for all valid points in the field */ struct stats cf; /* for composite fields, stats w/dup points ignored */ double has_cf; /* set if the cf stats are valid (unused) */ }; typedef struct stats *Stats; /* !!before you use this, read the comment following the data[1] line */ struct argblock { int len; /* length of this block including trailing bytes */ int isparallel; /* have we created a task group already? */ Field f; /* field to work on */ char *compname; /* component to compute stats for */ char *statname; /* stats component name */ char data[1]; /* compname and statname point here */ /* compname and statname point into the bytes directly following * this structure, which means you shouldn't try to to declare a * a variable of type struct argblock, only type struct argblock *. * (with a pointer, you can access the len and component name pointers * directly). all this silliness is necessary to run in parallel, * since all data for the tasks must be in global memory. * don't use sizeof() with this struct - use len instead. */ }; /* just too long for me to type, sorry */ #define ICH InvalidComponentHandle /* this should be in a system-wide header file (or it should be the * same number and i'm sure i use a different number in every file). */ #define MAXDIM 16 #if 0 static Field Compute_Wrapper(Field f, char *argblock, int argsize); static CompositeField CFCompute_Wrapper(CompositeField f, char *argblock, int argsize); #endif static Error supported(Array); static Stats Add__Stats(Object, char *, char *, int cfmember, struct stats *sp); static Error ComputeStats(Pointer); static Error AccumulateStats(Object, char *, char *, Stats, int *, int); static Error ConstantData(Array a, double *origin); static Error RegularData(Array a); static Error HasInvalid(Field f, char *, ICH *ih); static Error InvalidName(Field f, char *, char **); static Error buildblock(char *compname, char *statname, struct argblock **ab, int *len); /* temp */ static Error ProcessPartsXX(Object o, struct argblock *ab, int len); /* undocumented entry point in component.c */ extern Object Exists2(Object o, char *name, char *name2); /* the empty field call in libdx returns yes if there are components * but no positions component. this routine should be able to return * the data stats from fields without a position component. */ static int DXReallyEmptyField(Field f) { if (DXGetObjectClass((Object)f) != CLASS_FIELD) return 0; return DXGetEnumeratedComponentValue(f, 0, NULL) == NULL; } /* isnan doesn't exist for ibmpvs */ #if defined(ibmpvs) || defined(DXD_WIN) #define isnan(a) 0 #endif /* welcome to macro heaven. */ #define MAG2D(a) sqrt((double)(a)[0]*(a)[0] + (double)(a)[1]*(a)[1]) #define MAG3D(a) sqrt((double)(a)[0]*(a)[0] + \ (double)(a)[1]*(a)[1] + \ (double)(a)[2]*(a)[2]) #define DET2D(a) ((double)(a)[0]*(a)[3] - (double)(a)[1]*(a)[2]) #define DET3D(a) ((a)[0] * ((double)(a)[4]*(a)[8] - (double)(a)[5]*(a)[7]) - \ (a)[1] * ((double)(a)[3]*(a)[8] - (double)(a)[5]*(a)[6]) + \ (a)[2] * ((double)(a)[3]*(a)[7] - (double)(a)[4]*(a)[6])) static double determinant(double *a, int rank) { int i, j, k; int sign; double sum; double *minor, *dp; if (rank == 1) return a[0]; /* make space for the minor of the matrix (the matrix with the * i-th row and column removed). */ minor = (double *)DXAllocateLocal(sizeof(double) * (rank-1) * (rank-1)); if (!minor) return 0.0; for (i=0, sum=0.0, sign=1; imin = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = s->isvalid = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define vect2D_stats(type, tmin, tmax) \ static void \ type ## v2D_stat( type *dp, int size, struct stats *s, ICH inv) \ { \ int i; \ double thisval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define vect3D_stats(type, tmin, tmax) \ static void \ type ## v3D_stat( type *dp, int size, struct stats *s, ICH inv) \ { \ int i; \ double thisval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define vectND_stats(type, tmin, tmax) \ static void \ type ## vND_stat( type *dp, int size, int rank, struct stats *s, ICH inv) \ { \ int i, j; \ double thisval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define matr2D_stats(type, tmin, tmax) \ static void \ type ## m2D_stat( type *dp, int size, struct stats *s, ICH inv) \ { \ int i; \ double thisval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define matr3D_stats(type, tmin, tmax) \ static void \ type ## m3D_stat( type *dp, int size, struct stats *s, ICH inv) \ { \ int i; \ double thisval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define matrND_stats(type, tmin, tmax) \ static void \ type ## mND_stat( type *dp, int size, int rank, struct stats *s, ICH inv) \ { \ int i; \ double thisval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = 0; \ \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define regscalar_stats(type, tmin, tmax) \ static void \ type ## RS_stat( type *dp, type *inc, int size, struct stats *s, ICH inv) \ { \ int i; \ double thisval, startval, incval; \ \ s->min = tmax; s->max = tmin; s->sum = s->sumsq = 0.0; \ s->counts = s->isvalid = 0; \ \ startval = (double)*dp; \ incval = (double)*inc; \ for(i=0; i s->max) s->max = thisval; \ if(thisval < s->min) s->min = thisval; \ \ s->sum += thisval; \ s->sumsq += ((double)thisval * thisval); \ s->counts++; \ } \ \ if (s->counts > 0) { \ s->avg = s->sum / s->counts; \ s->isvalid = 1; \ } \ } #define scalar_conv(type) \ static void \ type ## S_conv( type *dp, int size, float *fp) \ { \ int i; \ for(i=0; iisvalid) goto error; sp.min = ap->min; sp.max = ap->max; sp.sum += ap->sum; sp.sumsq += ap->sumsq; sp.counts += ap->counts; } else { /* * if the component doesn't exist in any field, return error */ if (!DXExists (o, compname)) { DXSetError (ERROR_MISSING_DATA, "#10250", "input", compname); goto error; } #define STATLEN 12 /* strlen(" statistics"); */ statname = (char *)DXAllocateLocal (strlen(compname) + STATLEN); if (!statname) goto error; strcpy (statname, compname); strcat (statname, " statistics"); /* if the statistics don't already exist, compute them in parallel */ /* Exists2 is an undocumented call which requires that all fields * which contain the first component also contain the second. * when this section is done, all fields are guarenteed to have * a stats component. */ if (!Exists2 (o, compname, statname)) { /* * if we have to compute the stats, build a parm block * in global memory so we can use it with DXAddTask. * on a workstation we still need a place to store the * component name and the "compname statistics" strings, * so it isn't bad to do this even if this isn't an mp system. */ if (!buildblock(compname, statname, &ab, &len)) goto error; #if PP_TEST if ((_dxfProcessPartsG (o, CFCompute_Wrapper, ab, len, 0, 1)) == NULL) { if (DXGetError() != ERROR_NONE) goto error; /* else it's a valid empty field */ } #else /* if this is a simple field, compute the stats and return. * if this is a complicated object, schedule the stats to be * computed in parallel, and execute all the stats computations * with ExecuteTaskGroup. */ if ((ProcessPartsXX (o, ab, len) == ERROR)) goto error; if (ab->isparallel) { if (!DXExecuteTaskGroup()) goto error; } #endif #if 0 if ((DXProcessParts (o, Compute_Wrapper, ab, len, 0, 1)) == NULL) { if (DXGetError() != ERROR_NONE) goto error; /* else it's a valid empty field */ } #endif } /* when you get here, you know that all the fields which contain * a "compname" component also contain a "compname statistics" * component with valid stats. all that this routine needs to do * is combine all the stats components together to get the combined * min, max, avg, etc */ if (!AccumulateStats (o, compname, statname, &sp, &isvalid, 0)) goto error; /* if all the items in the input object are invalid, set this error * and return. it's unfortunate that it has to be an error, but * there is no way other than changing the interface to this routine * to indicate that the numbers coming back have nothing to do with * what the stats really are. */ if (!isvalid) { DXSetError (ERROR_INVALID_DATA, "#11905"); goto error; } } /* check for illegal values. since all values contribute to the sum * check it to see if it's NaN. if so, change all others to be also. */ if (isnan(sp.sum)) sp.min = sp.max = sp.sumsq = sp.avg = sp.sum; /* now that we know the values, fill in the ones the user asked for. * the temp variable is because on the pvs at one point we were having * round-off errors with very small numbers and the temp value was * getting slightly negative. asking for sqrt in that case segfaulted. */ if (min) *min = sp.min; if (max) *max = sp.max; if (avg && sp.counts) *avg = sp.sum / sp.counts; if (stdev && sp.counts > 1 && (sp.max != sp.min)) { double temp; temp = (sp.sumsq - (sp.sum * sp.sum / sp.counts)) / (sp.counts - 1); if (temp > 0.0) *stdev = sqrt(temp); } DXFree((Pointer)ab); DXFree((Pointer)statname); return OK; error: DXFree((Pointer)ab); DXFree((Pointer)statname); return ERROR; } /* * the library process parts routine doesn't stop at composite fields. * this is a replacement which does. */ static Error ProcessPartsXX(Object o, struct argblock *ab, int len) { Object subo, newo, subo2; Matrix m; Array a; Stats ap; int fixed, z; float pos; char *name; int i; switch (DXGetObjectClass(o)) { case CLASS_GROUP: if (ab->isparallel) { if (DXGetGroupClass((Group)o) == CLASS_COMPOSITEFIELD) { ab->f = (Field)o; if (!DXAddTask(ComputeStats, (Pointer)ab, len, 0.0)) goto error; } else { for (i=0; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) if (!ProcessPartsXX (subo, ab, len)) return ERROR; } } else { if (DXGetGroupClass((Group)o) == CLASS_COMPOSITEFIELD) { ab->f = (Field)o; return (ComputeStats((Pointer)ab)); } else { if (!DXCreateTaskGroup()) return ERROR; ab->isparallel = 1; for (i=0; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) if (!ProcessPartsXX (subo, ab, len)) return ERROR; } } break; case CLASS_FIELD: ab->f = (Field)o; if (ab->isparallel) { if (!DXAddTask(ComputeStats, (Pointer)ab, len, 0.0)) goto error; } else { if (!ComputeStats(ab)) return ERROR; } break; case CLASS_SCREEN: if (!DXGetScreenInfo((Screen)o, &subo, &fixed, &z)) return NULL; if (!ProcessPartsXX(subo, ab, len)) return ERROR; break; case CLASS_XFORM: if (!DXGetXformInfo((Xform)o, &subo, &m)) return NULL; if (!ProcessPartsXX(subo, ab, len)) return ERROR; break; case CLASS_CLIPPED: if (!DXGetClippedInfo((Clipped)o, &subo, &subo2)) return NULL; if (!ProcessPartsXX(subo, ab, len)) return ERROR; break; default: /* any other types of objects can't contain fields. */ break; } return OK; error: if (ab->isparallel) DXAbortTaskGroup(); return ERROR; } static Error buildblock(char *compname, char *statname, struct argblock **ab, int *len) { char *cp; *len = sizeof(struct argblock) + strlen(compname) + strlen(statname) + 2; *ab = (struct argblock *)DXAllocate(*len); if (!*ab) return ERROR; (*ab)->len = *len; (*ab)->isparallel = 0; (*ab)->f = NULL; cp = (*ab)->data; (*ab)->compname = cp; strcpy(cp, compname); cp += strlen(compname) + 1; (*ab)->statname = cp; strcpy(cp, statname); return OK; } /* * for ProcessParts - not used now. */ #if 0 static Field Compute_Wrapper(Field f, char *argblock, int argsize) { return ((Field)ComputeStats((Object)f, (struct argblock *)argblock, argsize)); } static CompositeField CFCompute_Wrapper(CompositeField f, char *argblock, int argsize) { return ((CompositeField)ComputeStats((Object)f, (struct argblock *)argblock, argsize)); } #endif /* * called only on fields or compositefields */ static Error ComputeStats(Pointer ptr) { Object o, subo; struct argblock *ab = (struct argblock *)ptr; int i; struct stats dummy_stats; /* should be only a field or composite field */ o = (Object)ab->f; if (!o) { DXSetError(ERROR_BAD_PARAMETER, "#11870"); return ERROR; } switch (DXGetObjectClass(o)) { case CLASS_GROUP: if (DXGetGroupClass((Group)o) != CLASS_COMPOSITEFIELD) goto error; /* partitioned data shares points, so compute it twice. once * ignoring duplicated points and once independently per field. * do the normal type of stats first, then do cf type. the order * is important here - there is some code at the end of Add__Stats * which cares (it wants the component to exist the next time around. */ #if PP_TEST if ((DXProcessParts (o, Compute_Wrapper, ab, len, 0, 1)) == NULL) return ERROR; #else for (i = 0; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) { if (DXReallyEmptyField((Field)subo)) continue; if (!Add__Stats(subo, ab->compname, ab->statname, 0,&dummy_stats)) return ERROR; } #endif /* this must be called on the entire composite field so it can * add the duplicated points to the invalid list. */ if (! DXInvalidateDupBoundary(o)) return ERROR; #if PP_TEST if ((DXProcessParts (o, Compute_Wrapper, ab, len, 0, 1)) == NULL) return ERROR; #else for (i = 0; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) { if (DXReallyEmptyField((Field)subo)) continue; if (!Add__Stats(subo, ab->compname, ab->statname, 1, &dummy_stats)) return ERROR; } #endif if (! DXRestoreDupBoundary(o)) return ERROR; break; case CLASS_FIELD: if (DXReallyEmptyField((Field)o)) break; if (!Add__Stats(o, ab->compname, ab->statname, 0, &dummy_stats)) return ERROR; break; default: error: DXSetError(ERROR_INTERNAL, "ComputeStats called with wrong object type"); return ERROR; } return OK; } static Error AccumulateStats(Object o, char *compname, char *statname, Stats sp, int *iv, int cfmember) { Object subo, newo, subo2; Matrix m; Array a; Stats ap; int fixed, z; float pos; char *name; int i; switch (DXGetObjectClass(o)) { case CLASS_GROUP: if (DXGetGroupClass((Group)o) == CLASS_COMPOSITEFIELD) { /* use the cf stats here */ for (i=0; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) if (!AccumulateStats(subo, compname, statname, sp, iv, 1)) return ERROR; } else { /* any other type of group uses normal stats */ for (i=0; subo = DXGetEnumeratedMember((Group)o, i, NULL); i++) if (!AccumulateStats(subo, compname, statname, sp, iv, 0)) return ERROR; } break; case CLASS_FIELD: /* ignore fields with no components. */ if (DXReallyEmptyField((Field)o)) break; a = (Array)DXGetComponentValue((Field)o, statname); if (!a) { a = (Array)DXGetComponentValue((Field)o, compname); if (!a) return OK; DXSetError(ERROR_MISSING_DATA, "#11880", compname); return ERROR; } ap = (struct stats *)DXGetArrayData(a); if (!ap) return ERROR; if (cfmember) ap++; if (ap->isvalid) { if (sp->min > ap->min) sp->min = ap->min; if (sp->max < ap->max) sp->max = ap->max; sp->sum += ap->sum; sp->sumsq += ap->sumsq; sp->counts += ap->counts; (*iv)++; } break; case CLASS_SCREEN: if (!DXGetScreenInfo((Screen)o, &subo, &fixed, &z)) return NULL; if (!AccumulateStats(subo, compname, statname, sp, iv, 0)) return ERROR; break; case CLASS_XFORM: if (!DXGetXformInfo((Xform)o, &subo, &m)) return NULL; if (!AccumulateStats(subo, compname, statname, sp, iv, 0)) return ERROR; break; case CLASS_CLIPPED: if (!DXGetClippedInfo((Clipped)o, &subo, &subo2)) return NULL; if (!AccumulateStats(subo, compname, statname, sp, iv, 0)) return ERROR; break; default: /* anything else can't have field objects below it, so it can * safely be passed through unchanged */ break; } return OK; } /* * On successful completion, fill in the Stats structure pointed to by sp. */ static Stats Add__Stats(Object o, char *compname, char *statname, int cfmember, struct stats *sptr) { Array a, na = NULL; Type type; Category category; int rank; Pointer dp, dp2; int isarray = 0; int i, j, size, rc = OK, shape[MAXDIM]; double origin[MAXDIM]; struct stats *s = NULL; ICH hasinvalid = NULL; /* this code either does a simple field, or a single array */ if (DXGetObjectClass(o) == CLASS_FIELD) { a = (Array)DXGetComponentValue((Field)o, compname); if (!a) { DXSetError(ERROR_MISSING_DATA, "#10250", "input", compname); return NULL; } /* are they already there? they must be if cfmember is set, they * may or may not be otherwise. */ na = (Array)DXGetComponentValue((Field)o, statname); if (cfmember && !na) { DXSetError(ERROR_INTERNAL, "problems computing partitioned statistics"); return NULL; } if (!cfmember && na) return (Stats)DXGetArrayData(na); na = NULL; if (!HasInvalid((Field)o, compname, &hasinvalid)) return NULL; } else { a = (Array)o; isarray++; } if (DXGetObjectClass((Object)a) != CLASS_ARRAY) { DXSetError(ERROR_BAD_CLASS, "#11890"); goto done1; } /* get size and shape */ DXGetArrayInfo(a, &size, &type, &category, &rank, shape); if (rank == 0) shape[0] = 1; /* check type, category, rank & shape */ if (supported(a) != OK) goto done1; s = sptr; memset(s, '\0', sizeof(struct stats)); /* after the stats array is set to 0, return if there is no data */ if (size <= 0) goto nodata; s->counts = size; /* if the data is constant, this converts the origins to double, and * does NOT call DXGetArrayData so the constant array doesn't get * expanded. */ if (ConstantData(a, origin) && !hasinvalid) { switch(shape[0]) { case 0: /* scalar == 1-vector == 1-matrix, etc */ case 1: doubleS_stat(origin, 1, s, NULL); break; case 2: /* 2-vector or 2x2 matrix */ switch(rank) { case 1: doublev2D_stat(origin, 1, s, NULL); break; case 2: doublem2D_stat(origin, 1, s, NULL); break; } break; case 3: /* 3-vector or 3x3 matrix */ switch(rank) { case 1: doublev3D_stat(origin, 1, s, NULL); break; case 2: doublem3D_stat(origin, 1, s, NULL); break; } break; default: /* N-vector or NxN matrix */ switch(rank) { case 1: doublevND_stat(origin, 1, shape[0], s, NULL); break; case 2: doublemND_stat(origin, 1, shape[0], s, NULL); break; } break; } /* fix sum and sum squared because we've only actually worked * on 1 data item instead of the real count. also have to * fix counts as well. */ s->sum = s->avg * size; s->sumsq = s->avg * s->avg * size; s->counts = size; if (s->counts <= 0) goto nodata; s->isvalid = 1; goto done; } /* define the macros used in the case statements below. */ #define EACHTYPE(ss) \ switch(type) { \ case TYPE_UBYTE: ubyte##ss##_stat((ubyte *)dp, size, s, hasinvalid); \ break; \ case TYPE_BYTE: byte##ss##_stat((byte *)dp, size, s, hasinvalid); \ break; \ case TYPE_USHORT: ushort##ss##_stat((ushort *)dp, size, s, hasinvalid); \ break; \ case TYPE_SHORT: short##ss##_stat((short *)dp, size, s, hasinvalid); \ break; \ case TYPE_UINT: uint##ss##_stat((uint *)dp, size, s, hasinvalid); \ break; \ case TYPE_INT: int##ss##_stat((int *)dp, size, s, hasinvalid); \ break; \ case TYPE_FLOAT: float##ss##_stat((float *)dp, size, s, hasinvalid); \ break; \ case TYPE_DOUBLE: double##ss##_stat((double *)dp, size, s, hasinvalid); \ break; \ } #define EACHTYPE_R(ss) \ switch(type) { \ case TYPE_UBYTE: \ ubyte##ss##_stat((ubyte *)dp, (ubyte *)dp2, size, s, hasinvalid); \ break; \ case TYPE_BYTE: \ byte##ss##_stat((byte *)dp, (byte *)dp2, size, s, hasinvalid); \ break; \ case TYPE_USHORT: \ ushort##ss##_stat((ushort *)dp, (ushort *)dp2, size, s, hasinvalid); \ break; \ case TYPE_SHORT: \ short##ss##_stat((short *)dp, (short *)dp2, size, s, hasinvalid); \ break; \ case TYPE_UINT: \ uint##ss##_stat((uint *)dp, (uint *)dp2, size, s, hasinvalid); \ break; \ case TYPE_INT: \ int##ss##_stat((int *)dp, (int *)dp2, size, s, hasinvalid); \ break; \ case TYPE_FLOAT: \ float##ss##_stat((float *)dp, (float *)dp2, size, s, hasinvalid); \ break; \ case TYPE_DOUBLE: \ double##ss##_stat((double *)dp, (double *)dp2, size, s, hasinvalid); \ break; \ } #define EACHTYPE_N(ss, len) \ switch(type) { \ case TYPE_UBYTE: \ ubyte##ss##_stat((ubyte *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_BYTE: \ byte##ss##_stat((byte *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_USHORT: \ ushort##ss##_stat((ushort *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_SHORT: \ short##ss##_stat((short *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_UINT: \ uint##ss##_stat((uint *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_INT: \ int##ss##_stat((int *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_FLOAT: \ float##ss##_stat((float *)dp, size, len, s, hasinvalid); \ break; \ case TYPE_DOUBLE: \ double##ss##_stat((double *)dp, size, len, s, hasinvalid); \ break; \ } /* if the data is scalar and regular, compute the stats without * calling DXGetArrayData so the data doesn't get expanded. */ if (RegularData(a)) { switch(shape[0]) { case 0: /* scalar == 1-vector == 1-matrix, etc */ case 1: dp = (Pointer)origin; dp2 = (Pointer)(origin+1); if (!DXGetRegularArrayInfo((RegularArray)a, NULL, dp, dp2)) goto error; EACHTYPE_R(RS); if (s->counts <= 0) goto nodata; s->isvalid = 1; goto done; default: /* for now, don't handle vectors or matricies */ break; } } /* must be irregular and/or have invalid parts */ dp = DXGetArrayData(a); if (!dp) goto error; switch(shape[0]) { case 0: /* scalar == 1-vector == 1-matrix, etc */ case 1: EACHTYPE(S); break; case 2: /* 2-vector or 2x2 matrix */ switch(rank) { case 1: EACHTYPE(v2D); break; case 2: EACHTYPE(m2D); break; } break; case 3: /* 3-vector or 3x3 matrix */ switch(rank) { case 1: EACHTYPE(v3D); break; case 2: EACHTYPE(m3D); break; } break; default: /* N-vector or NxN matrix */ switch(rank) { case 1: EACHTYPE_N(vND, shape[0]); break; case 2: EACHTYPE_N(mND, shape[0]); break; } break; } if (s->counts <= 0) goto nodata; #define STATLEN 12 /* strlen(" statistics"); */ nodata: /* fall thru */ done: if (!isarray) { int n = sizeof(struct stats) / sizeof(double); double one = 1.0; /* make new stats component if not composite field */ if (!cfmember) { na = DXNewArray(TYPE_DOUBLE, CATEGORY_REAL, 0); if(!na) goto error; if (!DXAddArrayData(na, 0, n, (Pointer)sptr)) goto error; /* if there was an invalid component, make the stats both der * on the data and on the invalid component. else just der data. */ if (hasinvalid) { char *cp; Array sl; if (!InvalidName((Field)o, compname, &cp)) goto error; sl = DXMakeStringList(2, compname, cp); DXFree(cp); if (!sl) goto error; if (!DXSetAttribute((Object)na, "der", (Object)sl)) { DXDelete((Object)sl); goto error; } } else if (!DXSetStringAttribute((Object)na, "der", compname)) goto error; if (!DXSetComponentValue((Field)o, statname, (Object)na)) goto error; } else { na = (Array)DXGetComponentValue((Field)o, statname); if (!na) goto error; if (!DXAddArrayData(na, n, n, (Pointer)sptr)) goto error; if (!DXAddArrayData(na, n*2, 1, &one)) goto error; } } done1: DXFreeInvalidComponentHandle(hasinvalid); return s; error: DXFreeInvalidComponentHandle(hasinvalid); DXDelete((Object)na); return NULL; } static Error supported(Array a) { Type type; Category category; int rank, shape[MAXDIM]; /* get size and shape */ DXGetArrayInfo(a, NULL, &type, &category, &rank, shape); /* recognized types */ switch(type) { case TYPE_UBYTE: case TYPE_BYTE: case TYPE_USHORT: case TYPE_SHORT: case TYPE_UINT: case TYPE_INT: case TYPE_FLOAT: case TYPE_DOUBLE: break; default: DXSetError(ERROR_NOT_IMPLEMENTED, "#11900", "type"); return ERROR; } /* scalar, all vectors, and square matricies */ if(rank > 2 || (rank == 2 && (shape[0] != shape[1]))) { DXSetError(ERROR_NOT_IMPLEMENTED, "#11900", "shape"); return ERROR; } if(category != CATEGORY_REAL) { DXSetError(ERROR_NOT_IMPLEMENTED, "#11900", "category"); return ERROR; } return OK; } Array DXScalarConvert(Array a) { int i, j, rc, nitems; int isregular = 0; int rank, shape[MAXDIM]; double origins[MAXDIM]; float norigin, ndelta; double x; Type type; Pointer dp, dp2; float *fp; Array na = NULL; /* type check: first object then datatype */ if (DXGetObjectClass((Object)a) != CLASS_ARRAY) { DXSetError(ERROR_INVALID_DATA, "#11920"); return NULL; } if (supported(a) != OK) return NULL; /* get the number of items in the input array */ if (!DXGetArrayInfo(a, &nitems, &type, NULL, &rank, shape)) return NULL; if (rank == 0) shape[0] = 1; if (nitems <= 0) { DXSetError(ERROR_INVALID_DATA, "#11910"); return NULL; } /* if the data is constant, create a rank 0, type float array * with the same number of items. the ConstantData() call * returns the data value as a double in the "origins" array, * so all the subsequent calls convert the double value to a float * with the XXX_conv() routines, and then finally norigin is used * when creating the scalar actual array. */ if (ConstantData(a, origins)) { switch(shape[0]) { case 0: /* scalars, 1-vectors and 1-matricies are same */ case 1: doubleS_conv(origins, 1, &norigin); break; case 2: /* 2-vectors or 2x2 matricies */ switch(rank) { case 1: doublev2D_conv(origins, 1, &norigin); break; case 2: doublem2D_conv(origins, 1, &norigin); break; } break; case 3: /* 3-vectors or 3x3 matricies */ switch(rank) { case 1: doublev3D_conv(origins, 1, &norigin); break; case 2: doublem3D_conv(origins, 1, &norigin); break; } break; default: /* N-vectors or NxN matricies */ switch(rank) { case 1: doublevND_conv(origins, 1, shape[0], &norigin); break; case 2: doublemND_conv(origins, 1, shape[0], &norigin); break; } break; } /* a constant array of scalars */ na = (Array)DXNewConstantArray(nitems, (Pointer)&norigin, TYPE_FLOAT, CATEGORY_REAL, 0); return na; } #undef EACHTYPE #define EACHTYPE(ss) \ switch(type) { \ case TYPE_UBYTE: ubyte##ss##_conv((ubyte *)dp, nitems, fp); \ break; \ case TYPE_BYTE: byte##ss##_conv((byte *)dp, nitems, fp); \ break; \ case TYPE_USHORT: ushort##ss##_conv((ushort *)dp, nitems, fp); \ break; \ case TYPE_SHORT: short##ss##_conv((short *)dp, nitems, fp); \ break; \ case TYPE_UINT: uint##ss##_conv((uint *)dp, nitems, fp); \ break; \ case TYPE_INT: int##ss##_conv((int *)dp, nitems, fp); \ break; \ case TYPE_FLOAT: float##ss##_conv((float *)dp, nitems, fp); \ break; \ case TYPE_DOUBLE: double##ss##_conv((double *)dp, nitems, fp); \ break; \ } #undef EACHTYPE_N #define EACHTYPE_N(ss, len) \ switch(type) { \ case TYPE_UBYTE: ubyte##ss##_conv((ubyte *)dp, nitems, len, fp); \ break; \ case TYPE_BYTE: byte##ss##_conv((byte *)dp, nitems, len, fp); \ break; \ case TYPE_USHORT: ushort##ss##_conv((ushort *)dp, nitems, len, fp); \ break; \ case TYPE_SHORT: short##ss##_conv((short *)dp, nitems, len, fp); \ break; \ case TYPE_UINT: uint##ss##_conv((uint *)dp, nitems, len, fp); \ break; \ case TYPE_INT: int##ss##_conv((int *)dp, nitems, len, fp); \ break; \ case TYPE_FLOAT: float##ss##_conv((float *)dp, nitems, len, fp); \ break; \ case TYPE_DOUBLE: double##ss##_conv((double *)dp, nitems, len, fp); \ break; \ } /* if the data is scalar and regular, compute the new float Regular array * without calling DXGetArrayData so the data doesn't get expanded. */ if (RegularData(a)) { switch(shape[0]) { case 0: /* scalar == 1-vector == 1-matrix, etc */ case 1: /* convert origin & deltas to float & create new reg array */ dp = (Pointer)origins; dp2 = (Pointer)(origins+1); if (!DXGetRegularArrayInfo((RegularArray)a, NULL, dp, dp2)) goto error; /* temporarily lie about the number of items */ i = nitems; nitems = 1; fp = (Pointer)&norigin; EACHTYPE(S); dp = (Pointer)(origins+1); fp = (Pointer)&ndelta; EACHTYPE(S); /* now restore them */ nitems = i; /* a regular array of scalars */ na = (Array)DXNewRegularArray(TYPE_FLOAT, 1, nitems, (Pointer)&norigin, (Pointer)&ndelta); return na; default: /* for now, don't handle vectors or matricies */ break; } } /* data must be irregular */ if (!(na = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0))) return NULL; if (!DXAddArrayData(na, nitems, 0, NULL)) goto error; dp = DXGetArrayData(a); if(!dp) goto error; fp = (float *)DXGetArrayData(na); if(!fp) goto error; switch(shape[0]) { case 0: /* scalars, 1-vectors and 1-matricies are same */ case 1: EACHTYPE(S); break; case 2: /* 2-vectors or 2x2 matricies */ switch(rank) { case 1: EACHTYPE(v2D); break; case 2: EACHTYPE(m2D); break; } break; case 3: /* 3-vectors or 3x3 matricies */ switch(rank) { case 1: EACHTYPE(v3D); break; case 2: EACHTYPE(m3D); break; } break; default: /* N-vectors or NxN matricies */ switch(rank) { case 1: EACHTYPE_N(vND, shape[0]); break; case 2: EACHTYPE_N(mND, shape[0]); break; } break; } /* return the new array */ return na; error: if (na) DXDelete((Object)na); return NULL; } /* returns OK if data is constant, and converts origins to double. * returns ERROR (but doesn't set an error code) if not constant data value. */ static Error ConstantData(Array a, double *origin) { Type type; Category cat; int rank; int i, shape[MAXDIM], count; ubyte *op; ubyte generic_origin[MAXDIM*sizeof(double)]; ubyte generic_delta[MAXDIM*sizeof(double)]; /* if scalar or vector, don't expand. for rank > 2, i'll need * to add more code. */ if (DXGetArrayClass(a) == CLASS_CONSTANTARRAY) { DXGetArrayInfo(a, NULL, &type, &cat, &rank, shape); if (cat != CATEGORY_REAL || rank >= 2) return ERROR; count = DXGetItemSize(a) / DXTypeSize(type); op = (ubyte *)DXGetConstantArrayData(a); #define MAKEDOUBLE(type) \ for (i=0; i= 2) return ERROR; if (rank == 1 && shape[0] != 1) return ERROR; return OK; } static Error HasInvalid(Field f, char *component, ICH *ih) { Error rc = OK; char *dep = NULL; char *invalid = NULL; *ih = NULL; if (!DXGetStringAttribute(DXGetComponentValue(f, component), "dep", &dep)) return OK; #define INVLEN 10 /* strlen("invalid "); */ if (!(invalid = (char *)DXAllocateLocal(strlen(dep) + INVLEN))) return ERROR; strcpy(invalid, "invalid "); strcat(invalid, dep); /* if no component, return NULL */ if (!DXGetComponentValue(f, invalid)) goto done; *ih = DXCreateInvalidComponentHandle((Object)f, NULL, dep); if (!*ih) { rc = ERROR; goto done; } done: DXFree((Pointer)invalid); return rc; } /* set the name of the invalid component which corresponds to the * given component. */ static Error InvalidName(Field f, char *component, char **invalid) { Error rc = OK; char *dep = NULL; *invalid = NULL; if (!DXGetStringAttribute(DXGetComponentValue(f, component), "dep", &dep)) return OK; #define INVLEN 10 /* strlen("invalid "); */ if (!(*invalid = (char *)DXAllocateLocal(strlen(dep) + INVLEN))) return ERROR; strcpy(*invalid, "invalid "); strcat(*invalid, dep); /* if no component, return ERROR - at the point we call this * routine, we know there are invalids. */ if (!DXGetComponentValue(f, *invalid)) { DXSetError(ERROR_INTERNAL, "looking for %s component, which is not found", *invalid); DXFree((Pointer)*invalid); *invalid = NULL; return ERROR; } return OK; }