Index: /trunk/psLib/src/dataIO/psFileUtilsErrors.dat =================================================================== --- /trunk/psLib/src/dataIO/psFileUtilsErrors.dat (revision 3270) +++ /trunk/psLib/src/dataIO/psFileUtilsErrors.dat (revision 3271) @@ -15,4 +15,5 @@ psLookupTable_DIVIDE_BY_ZERO Divide by zero error during interpolation. psLookupTable_INVALID_TYPE Invalid psLookupType, %d; +psLookupTable_TABLE_INVALID Lookup table is invalid. # psFits_NULL The input psFits object can not NULL. Index: /trunk/psLib/src/dataIO/psFileUtilsErrors.h =================================================================== --- /trunk/psLib/src/dataIO/psFileUtilsErrors.h (revision 3270) +++ /trunk/psLib/src/dataIO/psFileUtilsErrors.h (revision 3271) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-17 19:26:24 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 21:54:09 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -36,4 +36,5 @@ #define PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO "Divide by zero error during interpolation." #define PS_ERRORTEXT_psLookupTable_INVALID_TYPE "Invalid psLookupType, %d;" +#define PS_ERRORTEXT_psLookupTable_TABLE_INVALID "Lookup table is invalid." #define PS_ERRORTEXT_psFits_NULL "The input psFits object can not NULL." #define PS_ERRORTEXT_psFits_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" Index: /trunk/psLib/src/dataIO/psLookupTable.c =================================================================== --- /trunk/psLib/src/dataIO/psLookupTable.c (revision 3270) +++ /trunk/psLib/src/dataIO/psLookupTable.c (revision 3271) @@ -7,16 +7,21 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:24 $ +* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 21:54:09 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii */ #include #include #include +//#ifdef DARWIN #undef __STRICT_ANSI__ +//#endif #include -#include +//#ifdef DARWIN +#define __STRICT_ANSI__ +//#endif #include +#include #include "psMemory.h" @@ -438,4 +443,89 @@ } +#define UPDATE_VALID_TO_FROM(TABLE) \ +switch (TABLE->index->type.type) { \ +case PS_TYPE_U8: \ + TABLE->validFrom = (psF64)TABLE->index->data.U8[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U8[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S8: \ + TABLE->validFrom = (psF64)TABLE->index->data.S8[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S8[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U16: \ + TABLE->validFrom = (psF64)TABLE->index->data.U16[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U16[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S16: \ + TABLE->validFrom = (psF64)TABLE->index->data.S16[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S16[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U32: \ + TABLE->validFrom = (psF64)TABLE->index->data.U32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S32: \ + TABLE->validFrom = (psF64)TABLE->index->data.S32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U64: \ + TABLE->validFrom = (psF64)TABLE->index->data.U64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U64[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S64: \ + TABLE->validFrom = (psF64)TABLE->index->data.S64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S64[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_F32: \ + TABLE->validFrom = (psF64)TABLE->index->data.F32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.F32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_F64: \ + TABLE->validFrom = (psF64)TABLE->index->data.F64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.F64[TABLE->index->nalloc-1]; \ + break; \ +default: \ + TABLE->validFrom = (psF64)0; \ + TABLE->validTo = (psF64)0; \ + break; \ +} + +#define COPY_VECTOR_VALUES(VEC_OUT,INDEX_OUT,VEC_IN,INDEX_IN) \ +switch(((psVector*)(VEC_IN))->type.type) { \ +case PS_TYPE_U8: \ + ((psVector*)VEC_OUT)->data.U8[INDEX_OUT] = ((psVector*)VEC_IN)->data.U8[INDEX_IN]; \ + break; \ +case PS_TYPE_U16: \ + ((psVector*)VEC_OUT)->data.U16[INDEX_OUT] = ((psVector*)VEC_IN)->data.U16[INDEX_IN]; \ + break; \ +case PS_TYPE_U32: \ + ((psVector*)VEC_OUT)->data.U32[INDEX_OUT] = ((psVector*)VEC_IN)->data.U32[INDEX_IN]; \ + break; \ +case PS_TYPE_U64: \ + ((psVector*)VEC_OUT)->data.U64[INDEX_OUT] = ((psVector*)VEC_IN)->data.U64[INDEX_IN]; \ + break; \ +case PS_TYPE_S8: \ + ((psVector*)VEC_OUT)->data.S8[INDEX_OUT] = ((psVector*)VEC_IN)->data.S8[INDEX_IN]; \ + break; \ +case PS_TYPE_S16: \ + ((psVector*)VEC_OUT)->data.S16[INDEX_OUT] = ((psVector*)VEC_IN)->data.S16[INDEX_IN]; \ + break; \ +case PS_TYPE_S32: \ + ((psVector*)VEC_OUT)->data.S32[INDEX_OUT] = ((psVector*)VEC_IN)->data.S32[INDEX_IN]; \ + break; \ +case PS_TYPE_S64: \ + ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN]; \ + break; \ +case PS_TYPE_F32: \ + ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN]; \ + break; \ +case PS_TYPE_F64: \ + ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN]; \ + break; \ +default: \ + break; \ +} + + psLookupTable* psLookupTableRead(psLookupTable *table) { @@ -446,4 +536,5 @@ char *linePtr = NULL; psParseErrorType status = PS_PARSE_SUCCESS; + psParseErrorType lineStatus = PS_PARSE_SUCCESS; psU64 lineCount = 0; psU64 numRows = 0; @@ -451,18 +542,24 @@ psU64 failedLines = 0; FILE *fp = NULL; - psElemType elemType = PS_TYPE_PTR; + psElemType elemType; psVector *indexVec = NULL; psVector *valuesVec = NULL; psArray *values = NULL; - - - // Error checks + psBool sortIndex = false; + + // Check for NULL input table PS_PTR_CHECK_NULL(table,NULL); + + // Check for input table with NULL file name PS_PTR_CHECK_NULL(table->fileName,NULL); + + // Open table file specified by table->fileName if((fp=fopen(table->fileName, "r")) == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, table->fileName); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, + table->fileName); return table; } + // Initialize vector pointers indexVec = table->index; table->values = psArrayAlloc(1); @@ -473,5 +570,5 @@ // Loop through file to get numRows, numCols, and column data types - while(fgets(line, MAX_STRING_LENGTH, fp) != NULL) { + while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && (lineStatus == PS_PARSE_SUCCESS)) { // Initialize variables for new line @@ -485,5 +582,5 @@ // Determine column types from first line in data file after comments - while((strType=getToken(&linePtr, " ", &status)) != NULL) { + while((strType=getToken(&linePtr," ",&status))) { numCols++; typeLine = false; @@ -522,19 +619,21 @@ if(status) { - status = printError(lineCount, strValue, status); + printError(lineCount, strValue, status); failedLines++; + lineStatus = status; + status = PS_PARSE_SUCCESS; } psFree(strType); } } else { - // Parse and add values to all columns numRows++; numCols = 0; - while((strValue=getToken(&linePtr, " ", &status)) != NULL) { + while((strValue=getToken(&linePtr," ",&status))) { numCols++; // Realloc number of rows as you go if(numCols == 1) { + sortIndex = false; indexVec = psVectorRecycle(indexVec, numRows, indexVec->type.type); parseValue(indexVec, numRows-1, strValue, &status); @@ -546,6 +645,8 @@ if(status) { - status = printError(lineCount, strValue, status); + printError(lineCount, strValue, status); failedLines++; + lineStatus = status; + status = PS_PARSE_SUCCESS; } psFree(strValue); @@ -563,8 +664,192 @@ table->values = values; + // Set table values to indicate error detected during the read + if(lineStatus) { + table->numRows = 0; + table->numCols = 0; + table->validFrom = 0; + table->validTo = 0; + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psLookupTable_TABLE_INVALID); + } else { + // Check if index vector needs to be sorted + psVector* sortedIndex = psVectorAlloc(table->numRows,PS_TYPE_U32); + sortedIndex = psVectorSortIndex(sortedIndex,table->index); + for(psS32 i = 0; i < numRows; i++ ) { + if(sortedIndex->data.U32[i] != i) { + sortIndex = true; + break; + } + } + // Check if it is necessary to sort value vectors + if(sortIndex) { + // Allocate new index vector + psVector* newIndexVector = psVectorAlloc(numRows,indexVec->type.type); + // Allocate new value vectors + psArray* newValueArray = psArrayAlloc(numCols-1); + for(psS32 j = 0; j < numCols-1; j++) { + psS32 type = ((psVector*)(table->values->data[j]))->type.type; + newValueArray->data[j] = psVectorAlloc(numRows,type); + } + for(psS32 i = 0; i < numRows; i++) { + // Populate new index vector + psU32 sortIndex = sortedIndex->data.U32[i]; + COPY_VECTOR_VALUES(newIndexVector,i, indexVec,sortIndex) + // For every column populate new value vectors + for(psS32 j=0; j < numCols-1; j++) { + COPY_VECTOR_VALUES(newValueArray->data[j],i,table->values->data[j], sortIndex) + } + } + // Free old index vector + psFree(table->index); + // Free old value vectors + psFree(table->values); + // Assign new vector value to table + table->index = newIndexVector; + // Assign new value vectors to table array + table->values = newValueArray; + } + psFree(sortedIndex); + // Update the validTo and validFrom + UPDATE_VALID_TO_FROM(table) + } + fclose(fp); return table; } + +#define CONVERT_VALUE_TO_F64(VECTOR,INDEX,RESULT) \ +switch(VECTOR->type.type) { \ +case PS_TYPE_U8: \ + RESULT = (psF64)VECTOR->data.U8[INDEX]; \ + break; \ +case PS_TYPE_U16: \ + RESULT = (psF64)VECTOR->data.U16[INDEX]; \ + break; \ +case PS_TYPE_U32: \ + RESULT = (psF64)VECTOR->data.U32[INDEX]; \ + break; \ +case PS_TYPE_U64: \ + RESULT = (psF64)VECTOR->data.U64[INDEX]; \ + break; \ +case PS_TYPE_S8: \ + RESULT = (psF64)VECTOR->data.S8[INDEX]; \ + break; \ +case PS_TYPE_S16: \ + RESULT = (psF64)VECTOR->data.S16[INDEX]; \ + break; \ +case PS_TYPE_S32: \ + RESULT = (psF64)VECTOR->data.S32[INDEX]; \ + break; \ +case PS_TYPE_S64: \ + RESULT = (psF64)VECTOR->data.S64[INDEX]; \ + break; \ +case PS_TYPE_F32: \ + RESULT = (psF64)VECTOR->data.F32[INDEX]; \ + break; \ +case PS_TYPE_F64: \ + RESULT = VECTOR->data.F64[INDEX]; \ + break; \ +default: \ + RESULT = NAN; \ + break; \ +} + + +#define CHECK_LOWER_UPPER_BOUND(TABLE,INDEX,COLUMN) \ +switch (TABLE->index->type.type) { \ +case PS_TYPE_U8: \ + if( (psU8)index < TABLE->index->data.U8[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU8)index > TABLE->index->data.U8[numRows-1]) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U16: \ + if( (psU16)index < TABLE->index->data.U16[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU16)index > TABLE->index->data.U16[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U32: \ + if( (psU32)index < TABLE->index->data.U32[0]) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if ( (psU32)index > TABLE->index->data.U32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U64: \ + if( (psU64)index < TABLE->index->data.U64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU64)index > TABLE->index->data.U64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S8: \ + if( (psS8)index < TABLE->index->data.S8[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS8)index > TABLE->index->data.S8[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S16: \ + if( (psS16)index < TABLE->index->data.S16[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS16)index > TABLE->index->data.S16[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S32: \ + if( (psS32)index < TABLE->index->data.S32[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS32)index > TABLE->index->data.S32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S64: \ + if( (psS64)index < TABLE->index->data.S64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS64)index > TABLE->index->data.S64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_F32: \ + if( (psF32)index < TABLE->index->data.F32[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psF32)index > TABLE->index->data.F32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_F64: \ + if( index < TABLE->index->data.F64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( index > TABLE->index->data.F64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +default: \ + *status = PS_LOOKUP_ERROR; \ + return NAN; \ + break; \ +} \ +if(*status == PS_LOOKUP_PAST_TOP) { \ + CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),0,out) \ + return out; \ +} else if (*status == PS_LOOKUP_PAST_BOTTOM) { \ + CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),numRows-1,out) \ + return out; \ +} + psF64 psLookupTableInterpolate(psLookupTable *table, psF64 index, psU64 column, psLookupStatusType *status) @@ -576,63 +861,72 @@ psF64 out = 0.0; psF64 denom = 0.0; + psF64 convertVal = 0.0; + psF64 tempVal = 0.0; psVector *indexVec = NULL; psVector *valuesVec = NULL; psArray *values = NULL; - // Error checks - PS_PTR_CHECK_NULL(table,0.0); + // Set status to error prior to check since if checks fails it will immediately return + PS_PTR_CHECK_NULL(status,NAN); + *status = PS_LOOKUP_ERROR; + PS_PTR_CHECK_NULL(table,NAN); indexVec = table->index; values = table->values; numRows = table->numRows; numCols = table->numCols; - PS_PTR_CHECK_NULL(indexVec,0.0); - PS_PTR_CHECK_NULL(values,0.0); - PS_INT_CHECK_EQUALS(numRows, 0,0.0); - PS_INT_CHECK_EQUALS(numCols, 0,0.0); - PS_INT_CHECK_RANGE(column, 0, numCols-1, 0.0); - + PS_PTR_CHECK_NULL(indexVec,NAN); + PS_PTR_CHECK_NULL(values,NAN); + PS_INT_CHECK_EQUALS(numRows, 0,NAN); + PS_INT_CHECK_EQUALS(numCols, 0,NAN); + PS_INT_CHECK_RANGE(column, 0, (int)(numCols-1), NAN); valuesVec = (psVector*)values->data[column]; - PS_PTR_CHECK_NULL(indexVec,0.0); - - if(indexdata.F64[0] || indexvalidFrom) { - - // Index value past top of table - *status = PS_LOOKUP_PAST_TOP; - return valuesVec->data.F64[0]; - } else if(index>indexVec->data.F64[numRows-1] || index>table->validTo) { - - // Index value past bottom of table - *status = PS_LOOKUP_PAST_BOTTOM; - return valuesVec->data.F64[numRows-1]; + PS_PTR_CHECK_NULL(indexVec,NAN); + + // Set status to success since it passed all parameter checks + *status = PS_LOOKUP_SUCCESS; + + // Verify the index is within the bounds of the table + CHECK_LOWER_UPPER_BOUND(table,index,column) + + // Find location in table where specified index is between to entries + CONVERT_VALUE_TO_F64(indexVec, 0, convertVal) + while(index > convertVal ) { + hiIdx++; + if(hiIdx >= numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); + *status = PS_LOOKUP_ERROR; + return NAN; + } + CONVERT_VALUE_TO_F64(indexVec, hiIdx, convertVal) + } + + // Check for negative low index and generate error + loIdx = hiIdx--; + if(loIdx < 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); + *status = PS_LOOKUP_ERROR; + return NAN; + } + + // Perform linear interpolation to calculate return value + CONVERT_VALUE_TO_F64(indexVec, hiIdx, denom) + CONVERT_VALUE_TO_F64(indexVec, loIdx, convertVal); + denom -= convertVal; + if(fabs(denom) < FLT_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); + *status = PS_LOOKUP_ERROR; + return NAN; } else { - - // Interpolation - while(index > indexVec->data.F64[hiIdx]) { - hiIdx++; - if(hiIdx >= numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); - *status = PS_LOOKUP_ERROR; - return 0.0; - } - } - - loIdx = hiIdx--; - if(loIdx < 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); - *status = PS_LOOKUP_ERROR; - return 0.0; - } - - denom = indexVec->data.F64[hiIdx] - indexVec->data.F64[loIdx]; - if(fabs(denom) < FLT_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); - *status = PS_LOOKUP_ERROR; - return 0.0; - } else { - out = valuesVec->data.F64[loIdx]+(valuesVec->data.F64[hiIdx]-valuesVec->data.F64[loIdx]) - *(index-indexVec->data.F64[loIdx])/denom; - } + CONVERT_VALUE_TO_F64(valuesVec,hiIdx,tempVal) + CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) + tempVal -= convertVal; + CONVERT_VALUE_TO_F64(indexVec,loIdx,convertVal) + out = tempVal*(index-convertVal)/denom; + CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) + out += convertVal; } @@ -646,5 +940,4 @@ psVector *outVector = NULL; psLookupStatusType status = PS_LOOKUP_SUCCESS; - // Error checks @@ -654,5 +947,5 @@ PS_INT_CHECK_EQUALS(numCols, 0,NULL); - outVector = psVectorAlloc(numCols, PS_TYPE_F64); + outVector = psVectorAlloc(numCols+1, PS_TYPE_F64); // Fill vectors with results and status of results Index: /trunk/psLib/src/dataIO/psLookupTable.h =================================================================== --- /trunk/psLib/src/dataIO/psLookupTable.h (revision 3270) +++ /trunk/psLib/src/dataIO/psLookupTable.h (revision 3271) @@ -7,8 +7,8 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:24 $ +* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 21:54:09 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ Index: /trunk/psLib/src/fileUtils/psFileUtilsErrors.dat =================================================================== --- /trunk/psLib/src/fileUtils/psFileUtilsErrors.dat (revision 3270) +++ /trunk/psLib/src/fileUtils/psFileUtilsErrors.dat (revision 3271) @@ -15,4 +15,5 @@ psLookupTable_DIVIDE_BY_ZERO Divide by zero error during interpolation. psLookupTable_INVALID_TYPE Invalid psLookupType, %d; +psLookupTable_TABLE_INVALID Lookup table is invalid. # psFits_NULL The input psFits object can not NULL. Index: /trunk/psLib/src/fileUtils/psFileUtilsErrors.h =================================================================== --- /trunk/psLib/src/fileUtils/psFileUtilsErrors.h (revision 3270) +++ /trunk/psLib/src/fileUtils/psFileUtilsErrors.h (revision 3271) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-17 19:26:24 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 21:54:09 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -36,4 +36,5 @@ #define PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO "Divide by zero error during interpolation." #define PS_ERRORTEXT_psLookupTable_INVALID_TYPE "Invalid psLookupType, %d;" +#define PS_ERRORTEXT_psLookupTable_TABLE_INVALID "Lookup table is invalid." #define PS_ERRORTEXT_psFits_NULL "The input psFits object can not NULL." #define PS_ERRORTEXT_psFits_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" Index: /trunk/psLib/src/fileUtils/psLookupTable.c =================================================================== --- /trunk/psLib/src/fileUtils/psLookupTable.c (revision 3270) +++ /trunk/psLib/src/fileUtils/psLookupTable.c (revision 3271) @@ -7,16 +7,21 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:24 $ +* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 21:54:09 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii */ #include #include #include +//#ifdef DARWIN #undef __STRICT_ANSI__ +//#endif #include -#include +//#ifdef DARWIN +#define __STRICT_ANSI__ +//#endif #include +#include #include "psMemory.h" @@ -438,4 +443,89 @@ } +#define UPDATE_VALID_TO_FROM(TABLE) \ +switch (TABLE->index->type.type) { \ +case PS_TYPE_U8: \ + TABLE->validFrom = (psF64)TABLE->index->data.U8[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U8[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S8: \ + TABLE->validFrom = (psF64)TABLE->index->data.S8[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S8[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U16: \ + TABLE->validFrom = (psF64)TABLE->index->data.U16[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U16[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S16: \ + TABLE->validFrom = (psF64)TABLE->index->data.S16[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S16[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U32: \ + TABLE->validFrom = (psF64)TABLE->index->data.U32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S32: \ + TABLE->validFrom = (psF64)TABLE->index->data.S32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U64: \ + TABLE->validFrom = (psF64)TABLE->index->data.U64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U64[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S64: \ + TABLE->validFrom = (psF64)TABLE->index->data.S64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S64[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_F32: \ + TABLE->validFrom = (psF64)TABLE->index->data.F32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.F32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_F64: \ + TABLE->validFrom = (psF64)TABLE->index->data.F64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.F64[TABLE->index->nalloc-1]; \ + break; \ +default: \ + TABLE->validFrom = (psF64)0; \ + TABLE->validTo = (psF64)0; \ + break; \ +} + +#define COPY_VECTOR_VALUES(VEC_OUT,INDEX_OUT,VEC_IN,INDEX_IN) \ +switch(((psVector*)(VEC_IN))->type.type) { \ +case PS_TYPE_U8: \ + ((psVector*)VEC_OUT)->data.U8[INDEX_OUT] = ((psVector*)VEC_IN)->data.U8[INDEX_IN]; \ + break; \ +case PS_TYPE_U16: \ + ((psVector*)VEC_OUT)->data.U16[INDEX_OUT] = ((psVector*)VEC_IN)->data.U16[INDEX_IN]; \ + break; \ +case PS_TYPE_U32: \ + ((psVector*)VEC_OUT)->data.U32[INDEX_OUT] = ((psVector*)VEC_IN)->data.U32[INDEX_IN]; \ + break; \ +case PS_TYPE_U64: \ + ((psVector*)VEC_OUT)->data.U64[INDEX_OUT] = ((psVector*)VEC_IN)->data.U64[INDEX_IN]; \ + break; \ +case PS_TYPE_S8: \ + ((psVector*)VEC_OUT)->data.S8[INDEX_OUT] = ((psVector*)VEC_IN)->data.S8[INDEX_IN]; \ + break; \ +case PS_TYPE_S16: \ + ((psVector*)VEC_OUT)->data.S16[INDEX_OUT] = ((psVector*)VEC_IN)->data.S16[INDEX_IN]; \ + break; \ +case PS_TYPE_S32: \ + ((psVector*)VEC_OUT)->data.S32[INDEX_OUT] = ((psVector*)VEC_IN)->data.S32[INDEX_IN]; \ + break; \ +case PS_TYPE_S64: \ + ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN]; \ + break; \ +case PS_TYPE_F32: \ + ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN]; \ + break; \ +case PS_TYPE_F64: \ + ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN]; \ + break; \ +default: \ + break; \ +} + + psLookupTable* psLookupTableRead(psLookupTable *table) { @@ -446,4 +536,5 @@ char *linePtr = NULL; psParseErrorType status = PS_PARSE_SUCCESS; + psParseErrorType lineStatus = PS_PARSE_SUCCESS; psU64 lineCount = 0; psU64 numRows = 0; @@ -451,18 +542,24 @@ psU64 failedLines = 0; FILE *fp = NULL; - psElemType elemType = PS_TYPE_PTR; + psElemType elemType; psVector *indexVec = NULL; psVector *valuesVec = NULL; psArray *values = NULL; - - - // Error checks + psBool sortIndex = false; + + // Check for NULL input table PS_PTR_CHECK_NULL(table,NULL); + + // Check for input table with NULL file name PS_PTR_CHECK_NULL(table->fileName,NULL); + + // Open table file specified by table->fileName if((fp=fopen(table->fileName, "r")) == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, table->fileName); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, + table->fileName); return table; } + // Initialize vector pointers indexVec = table->index; table->values = psArrayAlloc(1); @@ -473,5 +570,5 @@ // Loop through file to get numRows, numCols, and column data types - while(fgets(line, MAX_STRING_LENGTH, fp) != NULL) { + while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && (lineStatus == PS_PARSE_SUCCESS)) { // Initialize variables for new line @@ -485,5 +582,5 @@ // Determine column types from first line in data file after comments - while((strType=getToken(&linePtr, " ", &status)) != NULL) { + while((strType=getToken(&linePtr," ",&status))) { numCols++; typeLine = false; @@ -522,19 +619,21 @@ if(status) { - status = printError(lineCount, strValue, status); + printError(lineCount, strValue, status); failedLines++; + lineStatus = status; + status = PS_PARSE_SUCCESS; } psFree(strType); } } else { - // Parse and add values to all columns numRows++; numCols = 0; - while((strValue=getToken(&linePtr, " ", &status)) != NULL) { + while((strValue=getToken(&linePtr," ",&status))) { numCols++; // Realloc number of rows as you go if(numCols == 1) { + sortIndex = false; indexVec = psVectorRecycle(indexVec, numRows, indexVec->type.type); parseValue(indexVec, numRows-1, strValue, &status); @@ -546,6 +645,8 @@ if(status) { - status = printError(lineCount, strValue, status); + printError(lineCount, strValue, status); failedLines++; + lineStatus = status; + status = PS_PARSE_SUCCESS; } psFree(strValue); @@ -563,8 +664,192 @@ table->values = values; + // Set table values to indicate error detected during the read + if(lineStatus) { + table->numRows = 0; + table->numCols = 0; + table->validFrom = 0; + table->validTo = 0; + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psLookupTable_TABLE_INVALID); + } else { + // Check if index vector needs to be sorted + psVector* sortedIndex = psVectorAlloc(table->numRows,PS_TYPE_U32); + sortedIndex = psVectorSortIndex(sortedIndex,table->index); + for(psS32 i = 0; i < numRows; i++ ) { + if(sortedIndex->data.U32[i] != i) { + sortIndex = true; + break; + } + } + // Check if it is necessary to sort value vectors + if(sortIndex) { + // Allocate new index vector + psVector* newIndexVector = psVectorAlloc(numRows,indexVec->type.type); + // Allocate new value vectors + psArray* newValueArray = psArrayAlloc(numCols-1); + for(psS32 j = 0; j < numCols-1; j++) { + psS32 type = ((psVector*)(table->values->data[j]))->type.type; + newValueArray->data[j] = psVectorAlloc(numRows,type); + } + for(psS32 i = 0; i < numRows; i++) { + // Populate new index vector + psU32 sortIndex = sortedIndex->data.U32[i]; + COPY_VECTOR_VALUES(newIndexVector,i, indexVec,sortIndex) + // For every column populate new value vectors + for(psS32 j=0; j < numCols-1; j++) { + COPY_VECTOR_VALUES(newValueArray->data[j],i,table->values->data[j], sortIndex) + } + } + // Free old index vector + psFree(table->index); + // Free old value vectors + psFree(table->values); + // Assign new vector value to table + table->index = newIndexVector; + // Assign new value vectors to table array + table->values = newValueArray; + } + psFree(sortedIndex); + // Update the validTo and validFrom + UPDATE_VALID_TO_FROM(table) + } + fclose(fp); return table; } + +#define CONVERT_VALUE_TO_F64(VECTOR,INDEX,RESULT) \ +switch(VECTOR->type.type) { \ +case PS_TYPE_U8: \ + RESULT = (psF64)VECTOR->data.U8[INDEX]; \ + break; \ +case PS_TYPE_U16: \ + RESULT = (psF64)VECTOR->data.U16[INDEX]; \ + break; \ +case PS_TYPE_U32: \ + RESULT = (psF64)VECTOR->data.U32[INDEX]; \ + break; \ +case PS_TYPE_U64: \ + RESULT = (psF64)VECTOR->data.U64[INDEX]; \ + break; \ +case PS_TYPE_S8: \ + RESULT = (psF64)VECTOR->data.S8[INDEX]; \ + break; \ +case PS_TYPE_S16: \ + RESULT = (psF64)VECTOR->data.S16[INDEX]; \ + break; \ +case PS_TYPE_S32: \ + RESULT = (psF64)VECTOR->data.S32[INDEX]; \ + break; \ +case PS_TYPE_S64: \ + RESULT = (psF64)VECTOR->data.S64[INDEX]; \ + break; \ +case PS_TYPE_F32: \ + RESULT = (psF64)VECTOR->data.F32[INDEX]; \ + break; \ +case PS_TYPE_F64: \ + RESULT = VECTOR->data.F64[INDEX]; \ + break; \ +default: \ + RESULT = NAN; \ + break; \ +} + + +#define CHECK_LOWER_UPPER_BOUND(TABLE,INDEX,COLUMN) \ +switch (TABLE->index->type.type) { \ +case PS_TYPE_U8: \ + if( (psU8)index < TABLE->index->data.U8[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU8)index > TABLE->index->data.U8[numRows-1]) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U16: \ + if( (psU16)index < TABLE->index->data.U16[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU16)index > TABLE->index->data.U16[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U32: \ + if( (psU32)index < TABLE->index->data.U32[0]) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if ( (psU32)index > TABLE->index->data.U32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U64: \ + if( (psU64)index < TABLE->index->data.U64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU64)index > TABLE->index->data.U64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S8: \ + if( (psS8)index < TABLE->index->data.S8[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS8)index > TABLE->index->data.S8[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S16: \ + if( (psS16)index < TABLE->index->data.S16[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS16)index > TABLE->index->data.S16[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S32: \ + if( (psS32)index < TABLE->index->data.S32[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS32)index > TABLE->index->data.S32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S64: \ + if( (psS64)index < TABLE->index->data.S64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS64)index > TABLE->index->data.S64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_F32: \ + if( (psF32)index < TABLE->index->data.F32[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psF32)index > TABLE->index->data.F32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_F64: \ + if( index < TABLE->index->data.F64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( index > TABLE->index->data.F64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +default: \ + *status = PS_LOOKUP_ERROR; \ + return NAN; \ + break; \ +} \ +if(*status == PS_LOOKUP_PAST_TOP) { \ + CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),0,out) \ + return out; \ +} else if (*status == PS_LOOKUP_PAST_BOTTOM) { \ + CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),numRows-1,out) \ + return out; \ +} + psF64 psLookupTableInterpolate(psLookupTable *table, psF64 index, psU64 column, psLookupStatusType *status) @@ -576,63 +861,72 @@ psF64 out = 0.0; psF64 denom = 0.0; + psF64 convertVal = 0.0; + psF64 tempVal = 0.0; psVector *indexVec = NULL; psVector *valuesVec = NULL; psArray *values = NULL; - // Error checks - PS_PTR_CHECK_NULL(table,0.0); + // Set status to error prior to check since if checks fails it will immediately return + PS_PTR_CHECK_NULL(status,NAN); + *status = PS_LOOKUP_ERROR; + PS_PTR_CHECK_NULL(table,NAN); indexVec = table->index; values = table->values; numRows = table->numRows; numCols = table->numCols; - PS_PTR_CHECK_NULL(indexVec,0.0); - PS_PTR_CHECK_NULL(values,0.0); - PS_INT_CHECK_EQUALS(numRows, 0,0.0); - PS_INT_CHECK_EQUALS(numCols, 0,0.0); - PS_INT_CHECK_RANGE(column, 0, numCols-1, 0.0); - + PS_PTR_CHECK_NULL(indexVec,NAN); + PS_PTR_CHECK_NULL(values,NAN); + PS_INT_CHECK_EQUALS(numRows, 0,NAN); + PS_INT_CHECK_EQUALS(numCols, 0,NAN); + PS_INT_CHECK_RANGE(column, 0, (int)(numCols-1), NAN); valuesVec = (psVector*)values->data[column]; - PS_PTR_CHECK_NULL(indexVec,0.0); - - if(indexdata.F64[0] || indexvalidFrom) { - - // Index value past top of table - *status = PS_LOOKUP_PAST_TOP; - return valuesVec->data.F64[0]; - } else if(index>indexVec->data.F64[numRows-1] || index>table->validTo) { - - // Index value past bottom of table - *status = PS_LOOKUP_PAST_BOTTOM; - return valuesVec->data.F64[numRows-1]; + PS_PTR_CHECK_NULL(indexVec,NAN); + + // Set status to success since it passed all parameter checks + *status = PS_LOOKUP_SUCCESS; + + // Verify the index is within the bounds of the table + CHECK_LOWER_UPPER_BOUND(table,index,column) + + // Find location in table where specified index is between to entries + CONVERT_VALUE_TO_F64(indexVec, 0, convertVal) + while(index > convertVal ) { + hiIdx++; + if(hiIdx >= numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); + *status = PS_LOOKUP_ERROR; + return NAN; + } + CONVERT_VALUE_TO_F64(indexVec, hiIdx, convertVal) + } + + // Check for negative low index and generate error + loIdx = hiIdx--; + if(loIdx < 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); + *status = PS_LOOKUP_ERROR; + return NAN; + } + + // Perform linear interpolation to calculate return value + CONVERT_VALUE_TO_F64(indexVec, hiIdx, denom) + CONVERT_VALUE_TO_F64(indexVec, loIdx, convertVal); + denom -= convertVal; + if(fabs(denom) < FLT_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); + *status = PS_LOOKUP_ERROR; + return NAN; } else { - - // Interpolation - while(index > indexVec->data.F64[hiIdx]) { - hiIdx++; - if(hiIdx >= numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); - *status = PS_LOOKUP_ERROR; - return 0.0; - } - } - - loIdx = hiIdx--; - if(loIdx < 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); - *status = PS_LOOKUP_ERROR; - return 0.0; - } - - denom = indexVec->data.F64[hiIdx] - indexVec->data.F64[loIdx]; - if(fabs(denom) < FLT_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); - *status = PS_LOOKUP_ERROR; - return 0.0; - } else { - out = valuesVec->data.F64[loIdx]+(valuesVec->data.F64[hiIdx]-valuesVec->data.F64[loIdx]) - *(index-indexVec->data.F64[loIdx])/denom; - } + CONVERT_VALUE_TO_F64(valuesVec,hiIdx,tempVal) + CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) + tempVal -= convertVal; + CONVERT_VALUE_TO_F64(indexVec,loIdx,convertVal) + out = tempVal*(index-convertVal)/denom; + CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) + out += convertVal; } @@ -646,5 +940,4 @@ psVector *outVector = NULL; psLookupStatusType status = PS_LOOKUP_SUCCESS; - // Error checks @@ -654,5 +947,5 @@ PS_INT_CHECK_EQUALS(numCols, 0,NULL); - outVector = psVectorAlloc(numCols, PS_TYPE_F64); + outVector = psVectorAlloc(numCols+1, PS_TYPE_F64); // Fill vectors with results and status of results Index: /trunk/psLib/src/fileUtils/psLookupTable.h =================================================================== --- /trunk/psLib/src/fileUtils/psLookupTable.h (revision 3270) +++ /trunk/psLib/src/fileUtils/psLookupTable.h (revision 3271) @@ -7,8 +7,8 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:24 $ +* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 21:54:09 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ Index: /trunk/psLib/src/types/psLookupTable.c =================================================================== --- /trunk/psLib/src/types/psLookupTable.c (revision 3270) +++ /trunk/psLib/src/types/psLookupTable.c (revision 3271) @@ -7,16 +7,21 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:24 $ +* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 21:54:09 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii */ #include #include #include +//#ifdef DARWIN #undef __STRICT_ANSI__ +//#endif #include -#include +//#ifdef DARWIN +#define __STRICT_ANSI__ +//#endif #include +#include #include "psMemory.h" @@ -438,4 +443,89 @@ } +#define UPDATE_VALID_TO_FROM(TABLE) \ +switch (TABLE->index->type.type) { \ +case PS_TYPE_U8: \ + TABLE->validFrom = (psF64)TABLE->index->data.U8[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U8[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S8: \ + TABLE->validFrom = (psF64)TABLE->index->data.S8[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S8[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U16: \ + TABLE->validFrom = (psF64)TABLE->index->data.U16[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U16[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S16: \ + TABLE->validFrom = (psF64)TABLE->index->data.S16[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S16[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U32: \ + TABLE->validFrom = (psF64)TABLE->index->data.U32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S32: \ + TABLE->validFrom = (psF64)TABLE->index->data.S32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_U64: \ + TABLE->validFrom = (psF64)TABLE->index->data.U64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.U64[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_S64: \ + TABLE->validFrom = (psF64)TABLE->index->data.S64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.S64[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_F32: \ + TABLE->validFrom = (psF64)TABLE->index->data.F32[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.F32[TABLE->index->nalloc-1]; \ + break; \ +case PS_TYPE_F64: \ + TABLE->validFrom = (psF64)TABLE->index->data.F64[0]; \ + TABLE->validTo = (psF64)TABLE->index->data.F64[TABLE->index->nalloc-1]; \ + break; \ +default: \ + TABLE->validFrom = (psF64)0; \ + TABLE->validTo = (psF64)0; \ + break; \ +} + +#define COPY_VECTOR_VALUES(VEC_OUT,INDEX_OUT,VEC_IN,INDEX_IN) \ +switch(((psVector*)(VEC_IN))->type.type) { \ +case PS_TYPE_U8: \ + ((psVector*)VEC_OUT)->data.U8[INDEX_OUT] = ((psVector*)VEC_IN)->data.U8[INDEX_IN]; \ + break; \ +case PS_TYPE_U16: \ + ((psVector*)VEC_OUT)->data.U16[INDEX_OUT] = ((psVector*)VEC_IN)->data.U16[INDEX_IN]; \ + break; \ +case PS_TYPE_U32: \ + ((psVector*)VEC_OUT)->data.U32[INDEX_OUT] = ((psVector*)VEC_IN)->data.U32[INDEX_IN]; \ + break; \ +case PS_TYPE_U64: \ + ((psVector*)VEC_OUT)->data.U64[INDEX_OUT] = ((psVector*)VEC_IN)->data.U64[INDEX_IN]; \ + break; \ +case PS_TYPE_S8: \ + ((psVector*)VEC_OUT)->data.S8[INDEX_OUT] = ((psVector*)VEC_IN)->data.S8[INDEX_IN]; \ + break; \ +case PS_TYPE_S16: \ + ((psVector*)VEC_OUT)->data.S16[INDEX_OUT] = ((psVector*)VEC_IN)->data.S16[INDEX_IN]; \ + break; \ +case PS_TYPE_S32: \ + ((psVector*)VEC_OUT)->data.S32[INDEX_OUT] = ((psVector*)VEC_IN)->data.S32[INDEX_IN]; \ + break; \ +case PS_TYPE_S64: \ + ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN]; \ + break; \ +case PS_TYPE_F32: \ + ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN]; \ + break; \ +case PS_TYPE_F64: \ + ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN]; \ + break; \ +default: \ + break; \ +} + + psLookupTable* psLookupTableRead(psLookupTable *table) { @@ -446,4 +536,5 @@ char *linePtr = NULL; psParseErrorType status = PS_PARSE_SUCCESS; + psParseErrorType lineStatus = PS_PARSE_SUCCESS; psU64 lineCount = 0; psU64 numRows = 0; @@ -451,18 +542,24 @@ psU64 failedLines = 0; FILE *fp = NULL; - psElemType elemType = PS_TYPE_PTR; + psElemType elemType; psVector *indexVec = NULL; psVector *valuesVec = NULL; psArray *values = NULL; - - - // Error checks + psBool sortIndex = false; + + // Check for NULL input table PS_PTR_CHECK_NULL(table,NULL); + + // Check for input table with NULL file name PS_PTR_CHECK_NULL(table->fileName,NULL); + + // Open table file specified by table->fileName if((fp=fopen(table->fileName, "r")) == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, table->fileName); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, + table->fileName); return table; } + // Initialize vector pointers indexVec = table->index; table->values = psArrayAlloc(1); @@ -473,5 +570,5 @@ // Loop through file to get numRows, numCols, and column data types - while(fgets(line, MAX_STRING_LENGTH, fp) != NULL) { + while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && (lineStatus == PS_PARSE_SUCCESS)) { // Initialize variables for new line @@ -485,5 +582,5 @@ // Determine column types from first line in data file after comments - while((strType=getToken(&linePtr, " ", &status)) != NULL) { + while((strType=getToken(&linePtr," ",&status))) { numCols++; typeLine = false; @@ -522,19 +619,21 @@ if(status) { - status = printError(lineCount, strValue, status); + printError(lineCount, strValue, status); failedLines++; + lineStatus = status; + status = PS_PARSE_SUCCESS; } psFree(strType); } } else { - // Parse and add values to all columns numRows++; numCols = 0; - while((strValue=getToken(&linePtr, " ", &status)) != NULL) { + while((strValue=getToken(&linePtr," ",&status))) { numCols++; // Realloc number of rows as you go if(numCols == 1) { + sortIndex = false; indexVec = psVectorRecycle(indexVec, numRows, indexVec->type.type); parseValue(indexVec, numRows-1, strValue, &status); @@ -546,6 +645,8 @@ if(status) { - status = printError(lineCount, strValue, status); + printError(lineCount, strValue, status); failedLines++; + lineStatus = status; + status = PS_PARSE_SUCCESS; } psFree(strValue); @@ -563,8 +664,192 @@ table->values = values; + // Set table values to indicate error detected during the read + if(lineStatus) { + table->numRows = 0; + table->numCols = 0; + table->validFrom = 0; + table->validTo = 0; + psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psLookupTable_TABLE_INVALID); + } else { + // Check if index vector needs to be sorted + psVector* sortedIndex = psVectorAlloc(table->numRows,PS_TYPE_U32); + sortedIndex = psVectorSortIndex(sortedIndex,table->index); + for(psS32 i = 0; i < numRows; i++ ) { + if(sortedIndex->data.U32[i] != i) { + sortIndex = true; + break; + } + } + // Check if it is necessary to sort value vectors + if(sortIndex) { + // Allocate new index vector + psVector* newIndexVector = psVectorAlloc(numRows,indexVec->type.type); + // Allocate new value vectors + psArray* newValueArray = psArrayAlloc(numCols-1); + for(psS32 j = 0; j < numCols-1; j++) { + psS32 type = ((psVector*)(table->values->data[j]))->type.type; + newValueArray->data[j] = psVectorAlloc(numRows,type); + } + for(psS32 i = 0; i < numRows; i++) { + // Populate new index vector + psU32 sortIndex = sortedIndex->data.U32[i]; + COPY_VECTOR_VALUES(newIndexVector,i, indexVec,sortIndex) + // For every column populate new value vectors + for(psS32 j=0; j < numCols-1; j++) { + COPY_VECTOR_VALUES(newValueArray->data[j],i,table->values->data[j], sortIndex) + } + } + // Free old index vector + psFree(table->index); + // Free old value vectors + psFree(table->values); + // Assign new vector value to table + table->index = newIndexVector; + // Assign new value vectors to table array + table->values = newValueArray; + } + psFree(sortedIndex); + // Update the validTo and validFrom + UPDATE_VALID_TO_FROM(table) + } + fclose(fp); return table; } + +#define CONVERT_VALUE_TO_F64(VECTOR,INDEX,RESULT) \ +switch(VECTOR->type.type) { \ +case PS_TYPE_U8: \ + RESULT = (psF64)VECTOR->data.U8[INDEX]; \ + break; \ +case PS_TYPE_U16: \ + RESULT = (psF64)VECTOR->data.U16[INDEX]; \ + break; \ +case PS_TYPE_U32: \ + RESULT = (psF64)VECTOR->data.U32[INDEX]; \ + break; \ +case PS_TYPE_U64: \ + RESULT = (psF64)VECTOR->data.U64[INDEX]; \ + break; \ +case PS_TYPE_S8: \ + RESULT = (psF64)VECTOR->data.S8[INDEX]; \ + break; \ +case PS_TYPE_S16: \ + RESULT = (psF64)VECTOR->data.S16[INDEX]; \ + break; \ +case PS_TYPE_S32: \ + RESULT = (psF64)VECTOR->data.S32[INDEX]; \ + break; \ +case PS_TYPE_S64: \ + RESULT = (psF64)VECTOR->data.S64[INDEX]; \ + break; \ +case PS_TYPE_F32: \ + RESULT = (psF64)VECTOR->data.F32[INDEX]; \ + break; \ +case PS_TYPE_F64: \ + RESULT = VECTOR->data.F64[INDEX]; \ + break; \ +default: \ + RESULT = NAN; \ + break; \ +} + + +#define CHECK_LOWER_UPPER_BOUND(TABLE,INDEX,COLUMN) \ +switch (TABLE->index->type.type) { \ +case PS_TYPE_U8: \ + if( (psU8)index < TABLE->index->data.U8[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU8)index > TABLE->index->data.U8[numRows-1]) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U16: \ + if( (psU16)index < TABLE->index->data.U16[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU16)index > TABLE->index->data.U16[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U32: \ + if( (psU32)index < TABLE->index->data.U32[0]) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if ( (psU32)index > TABLE->index->data.U32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_U64: \ + if( (psU64)index < TABLE->index->data.U64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psU64)index > TABLE->index->data.U64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S8: \ + if( (psS8)index < TABLE->index->data.S8[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS8)index > TABLE->index->data.S8[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S16: \ + if( (psS16)index < TABLE->index->data.S16[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS16)index > TABLE->index->data.S16[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S32: \ + if( (psS32)index < TABLE->index->data.S32[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS32)index > TABLE->index->data.S32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_S64: \ + if( (psS64)index < TABLE->index->data.S64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psS64)index > TABLE->index->data.S64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_F32: \ + if( (psF32)index < TABLE->index->data.F32[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( (psF32)index > TABLE->index->data.F32[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +case PS_TYPE_F64: \ + if( index < TABLE->index->data.F64[0] ) { \ + *status = PS_LOOKUP_PAST_TOP; \ + } \ + if( index > TABLE->index->data.F64[numRows-1] ) { \ + *status = PS_LOOKUP_PAST_BOTTOM; \ + } \ + break; \ +default: \ + *status = PS_LOOKUP_ERROR; \ + return NAN; \ + break; \ +} \ +if(*status == PS_LOOKUP_PAST_TOP) { \ + CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),0,out) \ + return out; \ +} else if (*status == PS_LOOKUP_PAST_BOTTOM) { \ + CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),numRows-1,out) \ + return out; \ +} + psF64 psLookupTableInterpolate(psLookupTable *table, psF64 index, psU64 column, psLookupStatusType *status) @@ -576,63 +861,72 @@ psF64 out = 0.0; psF64 denom = 0.0; + psF64 convertVal = 0.0; + psF64 tempVal = 0.0; psVector *indexVec = NULL; psVector *valuesVec = NULL; psArray *values = NULL; - // Error checks - PS_PTR_CHECK_NULL(table,0.0); + // Set status to error prior to check since if checks fails it will immediately return + PS_PTR_CHECK_NULL(status,NAN); + *status = PS_LOOKUP_ERROR; + PS_PTR_CHECK_NULL(table,NAN); indexVec = table->index; values = table->values; numRows = table->numRows; numCols = table->numCols; - PS_PTR_CHECK_NULL(indexVec,0.0); - PS_PTR_CHECK_NULL(values,0.0); - PS_INT_CHECK_EQUALS(numRows, 0,0.0); - PS_INT_CHECK_EQUALS(numCols, 0,0.0); - PS_INT_CHECK_RANGE(column, 0, numCols-1, 0.0); - + PS_PTR_CHECK_NULL(indexVec,NAN); + PS_PTR_CHECK_NULL(values,NAN); + PS_INT_CHECK_EQUALS(numRows, 0,NAN); + PS_INT_CHECK_EQUALS(numCols, 0,NAN); + PS_INT_CHECK_RANGE(column, 0, (int)(numCols-1), NAN); valuesVec = (psVector*)values->data[column]; - PS_PTR_CHECK_NULL(indexVec,0.0); - - if(indexdata.F64[0] || indexvalidFrom) { - - // Index value past top of table - *status = PS_LOOKUP_PAST_TOP; - return valuesVec->data.F64[0]; - } else if(index>indexVec->data.F64[numRows-1] || index>table->validTo) { - - // Index value past bottom of table - *status = PS_LOOKUP_PAST_BOTTOM; - return valuesVec->data.F64[numRows-1]; + PS_PTR_CHECK_NULL(indexVec,NAN); + + // Set status to success since it passed all parameter checks + *status = PS_LOOKUP_SUCCESS; + + // Verify the index is within the bounds of the table + CHECK_LOWER_UPPER_BOUND(table,index,column) + + // Find location in table where specified index is between to entries + CONVERT_VALUE_TO_F64(indexVec, 0, convertVal) + while(index > convertVal ) { + hiIdx++; + if(hiIdx >= numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); + *status = PS_LOOKUP_ERROR; + return NAN; + } + CONVERT_VALUE_TO_F64(indexVec, hiIdx, convertVal) + } + + // Check for negative low index and generate error + loIdx = hiIdx--; + if(loIdx < 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); + *status = PS_LOOKUP_ERROR; + return NAN; + } + + // Perform linear interpolation to calculate return value + CONVERT_VALUE_TO_F64(indexVec, hiIdx, denom) + CONVERT_VALUE_TO_F64(indexVec, loIdx, convertVal); + denom -= convertVal; + if(fabs(denom) < FLT_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); + *status = PS_LOOKUP_ERROR; + return NAN; } else { - - // Interpolation - while(index > indexVec->data.F64[hiIdx]) { - hiIdx++; - if(hiIdx >= numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); - *status = PS_LOOKUP_ERROR; - return 0.0; - } - } - - loIdx = hiIdx--; - if(loIdx < 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); - *status = PS_LOOKUP_ERROR; - return 0.0; - } - - denom = indexVec->data.F64[hiIdx] - indexVec->data.F64[loIdx]; - if(fabs(denom) < FLT_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); - *status = PS_LOOKUP_ERROR; - return 0.0; - } else { - out = valuesVec->data.F64[loIdx]+(valuesVec->data.F64[hiIdx]-valuesVec->data.F64[loIdx]) - *(index-indexVec->data.F64[loIdx])/denom; - } + CONVERT_VALUE_TO_F64(valuesVec,hiIdx,tempVal) + CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) + tempVal -= convertVal; + CONVERT_VALUE_TO_F64(indexVec,loIdx,convertVal) + out = tempVal*(index-convertVal)/denom; + CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) + out += convertVal; } @@ -646,5 +940,4 @@ psVector *outVector = NULL; psLookupStatusType status = PS_LOOKUP_SUCCESS; - // Error checks @@ -654,5 +947,5 @@ PS_INT_CHECK_EQUALS(numCols, 0,NULL); - outVector = psVectorAlloc(numCols, PS_TYPE_F64); + outVector = psVectorAlloc(numCols+1, PS_TYPE_F64); // Fill vectors with results and status of results Index: /trunk/psLib/src/types/psLookupTable.h =================================================================== --- /trunk/psLib/src/types/psLookupTable.h (revision 3270) +++ /trunk/psLib/src/types/psLookupTable.h (revision 3271) @@ -7,8 +7,8 @@ * @author Ross Harman, MHPCC * -* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-02-17 19:26:24 $ +* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 21:54:09 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004 Maui High Performance Computing Center, University of Hawaii */