Changeset 3271
- Timestamp:
- Feb 17, 2005, 11:54:09 AM (21 years ago)
- Location:
- trunk/psLib/src
- Files:
-
- 10 edited
-
dataIO/psFileUtilsErrors.dat (modified) (1 diff)
-
dataIO/psFileUtilsErrors.h (modified) (2 diffs)
-
dataIO/psLookupTable.c (modified) (12 diffs)
-
dataIO/psLookupTable.h (modified) (1 diff)
-
fileUtils/psFileUtilsErrors.dat (modified) (1 diff)
-
fileUtils/psFileUtilsErrors.h (modified) (2 diffs)
-
fileUtils/psLookupTable.c (modified) (12 diffs)
-
fileUtils/psLookupTable.h (modified) (1 diff)
-
types/psLookupTable.c (modified) (12 diffs)
-
types/psLookupTable.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataIO/psFileUtilsErrors.dat
r3028 r3271 15 15 psLookupTable_DIVIDE_BY_ZERO Divide by zero error during interpolation. 16 16 psLookupTable_INVALID_TYPE Invalid psLookupType, %d; 17 psLookupTable_TABLE_INVALID Lookup table is invalid. 17 18 # 18 19 psFits_NULL The input psFits object can not NULL. -
trunk/psLib/src/dataIO/psFileUtilsErrors.h
r3264 r3271 7 7 * @author Robert DeSonia, MHPCC 8 8 * 9 * @version $Revision: 1.1 1$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 36 36 #define PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO "Divide by zero error during interpolation." 37 37 #define PS_ERRORTEXT_psLookupTable_INVALID_TYPE "Invalid psLookupType, %d;" 38 #define PS_ERRORTEXT_psLookupTable_TABLE_INVALID "Lookup table is invalid." 38 39 #define PS_ERRORTEXT_psFits_NULL "The input psFits object can not NULL." 39 40 #define PS_ERRORTEXT_psFits_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" -
trunk/psLib/src/dataIO/psLookupTable.c
r3264 r3271 7 7 * @author Ross Harman, MHPCC 8 8 * 9 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 * Copyright 2004- 2005 Maui High Performance Computing Center, University of Hawaii12 * Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii 13 13 */ 14 14 #include <stdio.h> 15 15 #include <string.h> 16 16 #include <ctype.h> 17 //#ifdef DARWIN 17 18 #undef __STRICT_ANSI__ 19 //#endif 18 20 #include <stdlib.h> 19 #include <limits.h> 21 //#ifdef DARWIN 22 #define __STRICT_ANSI__ 23 //#endif 20 24 #include <math.h> 25 #include <stdlib.h> 21 26 22 27 #include "psMemory.h" … … 438 443 } 439 444 445 #define UPDATE_VALID_TO_FROM(TABLE) \ 446 switch (TABLE->index->type.type) { \ 447 case PS_TYPE_U8: \ 448 TABLE->validFrom = (psF64)TABLE->index->data.U8[0]; \ 449 TABLE->validTo = (psF64)TABLE->index->data.U8[TABLE->index->nalloc-1]; \ 450 break; \ 451 case PS_TYPE_S8: \ 452 TABLE->validFrom = (psF64)TABLE->index->data.S8[0]; \ 453 TABLE->validTo = (psF64)TABLE->index->data.S8[TABLE->index->nalloc-1]; \ 454 break; \ 455 case PS_TYPE_U16: \ 456 TABLE->validFrom = (psF64)TABLE->index->data.U16[0]; \ 457 TABLE->validTo = (psF64)TABLE->index->data.U16[TABLE->index->nalloc-1]; \ 458 break; \ 459 case PS_TYPE_S16: \ 460 TABLE->validFrom = (psF64)TABLE->index->data.S16[0]; \ 461 TABLE->validTo = (psF64)TABLE->index->data.S16[TABLE->index->nalloc-1]; \ 462 break; \ 463 case PS_TYPE_U32: \ 464 TABLE->validFrom = (psF64)TABLE->index->data.U32[0]; \ 465 TABLE->validTo = (psF64)TABLE->index->data.U32[TABLE->index->nalloc-1]; \ 466 break; \ 467 case PS_TYPE_S32: \ 468 TABLE->validFrom = (psF64)TABLE->index->data.S32[0]; \ 469 TABLE->validTo = (psF64)TABLE->index->data.S32[TABLE->index->nalloc-1]; \ 470 break; \ 471 case PS_TYPE_U64: \ 472 TABLE->validFrom = (psF64)TABLE->index->data.U64[0]; \ 473 TABLE->validTo = (psF64)TABLE->index->data.U64[TABLE->index->nalloc-1]; \ 474 break; \ 475 case PS_TYPE_S64: \ 476 TABLE->validFrom = (psF64)TABLE->index->data.S64[0]; \ 477 TABLE->validTo = (psF64)TABLE->index->data.S64[TABLE->index->nalloc-1]; \ 478 break; \ 479 case PS_TYPE_F32: \ 480 TABLE->validFrom = (psF64)TABLE->index->data.F32[0]; \ 481 TABLE->validTo = (psF64)TABLE->index->data.F32[TABLE->index->nalloc-1]; \ 482 break; \ 483 case PS_TYPE_F64: \ 484 TABLE->validFrom = (psF64)TABLE->index->data.F64[0]; \ 485 TABLE->validTo = (psF64)TABLE->index->data.F64[TABLE->index->nalloc-1]; \ 486 break; \ 487 default: \ 488 TABLE->validFrom = (psF64)0; \ 489 TABLE->validTo = (psF64)0; \ 490 break; \ 491 } 492 493 #define COPY_VECTOR_VALUES(VEC_OUT,INDEX_OUT,VEC_IN,INDEX_IN) \ 494 switch(((psVector*)(VEC_IN))->type.type) { \ 495 case PS_TYPE_U8: \ 496 ((psVector*)VEC_OUT)->data.U8[INDEX_OUT] = ((psVector*)VEC_IN)->data.U8[INDEX_IN]; \ 497 break; \ 498 case PS_TYPE_U16: \ 499 ((psVector*)VEC_OUT)->data.U16[INDEX_OUT] = ((psVector*)VEC_IN)->data.U16[INDEX_IN]; \ 500 break; \ 501 case PS_TYPE_U32: \ 502 ((psVector*)VEC_OUT)->data.U32[INDEX_OUT] = ((psVector*)VEC_IN)->data.U32[INDEX_IN]; \ 503 break; \ 504 case PS_TYPE_U64: \ 505 ((psVector*)VEC_OUT)->data.U64[INDEX_OUT] = ((psVector*)VEC_IN)->data.U64[INDEX_IN]; \ 506 break; \ 507 case PS_TYPE_S8: \ 508 ((psVector*)VEC_OUT)->data.S8[INDEX_OUT] = ((psVector*)VEC_IN)->data.S8[INDEX_IN]; \ 509 break; \ 510 case PS_TYPE_S16: \ 511 ((psVector*)VEC_OUT)->data.S16[INDEX_OUT] = ((psVector*)VEC_IN)->data.S16[INDEX_IN]; \ 512 break; \ 513 case PS_TYPE_S32: \ 514 ((psVector*)VEC_OUT)->data.S32[INDEX_OUT] = ((psVector*)VEC_IN)->data.S32[INDEX_IN]; \ 515 break; \ 516 case PS_TYPE_S64: \ 517 ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN]; \ 518 break; \ 519 case PS_TYPE_F32: \ 520 ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN]; \ 521 break; \ 522 case PS_TYPE_F64: \ 523 ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN]; \ 524 break; \ 525 default: \ 526 break; \ 527 } 528 529 440 530 psLookupTable* psLookupTableRead(psLookupTable *table) 441 531 { … … 446 536 char *linePtr = NULL; 447 537 psParseErrorType status = PS_PARSE_SUCCESS; 538 psParseErrorType lineStatus = PS_PARSE_SUCCESS; 448 539 psU64 lineCount = 0; 449 540 psU64 numRows = 0; … … 451 542 psU64 failedLines = 0; 452 543 FILE *fp = NULL; 453 psElemType elemType = PS_TYPE_PTR;544 psElemType elemType; 454 545 psVector *indexVec = NULL; 455 546 psVector *valuesVec = NULL; 456 547 psArray *values = NULL; 457 458 459 // Error checks548 psBool sortIndex = false; 549 550 // Check for NULL input table 460 551 PS_PTR_CHECK_NULL(table,NULL); 552 553 // Check for input table with NULL file name 461 554 PS_PTR_CHECK_NULL(table->fileName,NULL); 555 556 // Open table file specified by table->fileName 462 557 if((fp=fopen(table->fileName, "r")) == NULL) { 463 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, table->fileName); 558 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, 559 table->fileName); 464 560 return table; 465 561 } 466 562 563 // Initialize vector pointers 467 564 indexVec = table->index; 468 565 table->values = psArrayAlloc(1); … … 473 570 474 571 // Loop through file to get numRows, numCols, and column data types 475 while( fgets(line, MAX_STRING_LENGTH, fp) != NULL) {572 while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && (lineStatus == PS_PARSE_SUCCESS)) { 476 573 477 574 // Initialize variables for new line … … 485 582 486 583 // Determine column types from first line in data file after comments 487 while((strType=getToken(&linePtr, " ", &status)) != NULL) {584 while((strType=getToken(&linePtr," ",&status))) { 488 585 numCols++; 489 586 typeLine = false; … … 522 619 523 620 if(status) { 524 status =printError(lineCount, strValue, status);621 printError(lineCount, strValue, status); 525 622 failedLines++; 623 lineStatus = status; 624 status = PS_PARSE_SUCCESS; 526 625 } 527 626 psFree(strType); 528 627 } 529 628 } else { 530 531 629 // Parse and add values to all columns 532 630 numRows++; 533 631 numCols = 0; 534 while((strValue=getToken(&linePtr, " ", &status)) != NULL) {632 while((strValue=getToken(&linePtr," ",&status))) { 535 633 numCols++; 536 634 537 635 // Realloc number of rows as you go 538 636 if(numCols == 1) { 637 sortIndex = false; 539 638 indexVec = psVectorRecycle(indexVec, numRows, indexVec->type.type); 540 639 parseValue(indexVec, numRows-1, strValue, &status); … … 546 645 547 646 if(status) { 548 status =printError(lineCount, strValue, status);647 printError(lineCount, strValue, status); 549 648 failedLines++; 649 lineStatus = status; 650 status = PS_PARSE_SUCCESS; 550 651 } 551 652 psFree(strValue); … … 563 664 table->values = values; 564 665 666 // Set table values to indicate error detected during the read 667 if(lineStatus) { 668 table->numRows = 0; 669 table->numCols = 0; 670 table->validFrom = 0; 671 table->validTo = 0; 672 psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psLookupTable_TABLE_INVALID); 673 } else { 674 // Check if index vector needs to be sorted 675 psVector* sortedIndex = psVectorAlloc(table->numRows,PS_TYPE_U32); 676 sortedIndex = psVectorSortIndex(sortedIndex,table->index); 677 for(psS32 i = 0; i < numRows; i++ ) { 678 if(sortedIndex->data.U32[i] != i) { 679 sortIndex = true; 680 break; 681 } 682 } 683 // Check if it is necessary to sort value vectors 684 if(sortIndex) { 685 // Allocate new index vector 686 psVector* newIndexVector = psVectorAlloc(numRows,indexVec->type.type); 687 // Allocate new value vectors 688 psArray* newValueArray = psArrayAlloc(numCols-1); 689 for(psS32 j = 0; j < numCols-1; j++) { 690 psS32 type = ((psVector*)(table->values->data[j]))->type.type; 691 newValueArray->data[j] = psVectorAlloc(numRows,type); 692 } 693 for(psS32 i = 0; i < numRows; i++) { 694 // Populate new index vector 695 psU32 sortIndex = sortedIndex->data.U32[i]; 696 COPY_VECTOR_VALUES(newIndexVector,i, indexVec,sortIndex) 697 // For every column populate new value vectors 698 for(psS32 j=0; j < numCols-1; j++) { 699 COPY_VECTOR_VALUES(newValueArray->data[j],i,table->values->data[j], sortIndex) 700 } 701 } 702 // Free old index vector 703 psFree(table->index); 704 // Free old value vectors 705 psFree(table->values); 706 // Assign new vector value to table 707 table->index = newIndexVector; 708 // Assign new value vectors to table array 709 table->values = newValueArray; 710 } 711 psFree(sortedIndex); 712 // Update the validTo and validFrom 713 UPDATE_VALID_TO_FROM(table) 714 } 715 565 716 fclose(fp); 566 717 567 718 return table; 568 719 } 720 721 #define CONVERT_VALUE_TO_F64(VECTOR,INDEX,RESULT) \ 722 switch(VECTOR->type.type) { \ 723 case PS_TYPE_U8: \ 724 RESULT = (psF64)VECTOR->data.U8[INDEX]; \ 725 break; \ 726 case PS_TYPE_U16: \ 727 RESULT = (psF64)VECTOR->data.U16[INDEX]; \ 728 break; \ 729 case PS_TYPE_U32: \ 730 RESULT = (psF64)VECTOR->data.U32[INDEX]; \ 731 break; \ 732 case PS_TYPE_U64: \ 733 RESULT = (psF64)VECTOR->data.U64[INDEX]; \ 734 break; \ 735 case PS_TYPE_S8: \ 736 RESULT = (psF64)VECTOR->data.S8[INDEX]; \ 737 break; \ 738 case PS_TYPE_S16: \ 739 RESULT = (psF64)VECTOR->data.S16[INDEX]; \ 740 break; \ 741 case PS_TYPE_S32: \ 742 RESULT = (psF64)VECTOR->data.S32[INDEX]; \ 743 break; \ 744 case PS_TYPE_S64: \ 745 RESULT = (psF64)VECTOR->data.S64[INDEX]; \ 746 break; \ 747 case PS_TYPE_F32: \ 748 RESULT = (psF64)VECTOR->data.F32[INDEX]; \ 749 break; \ 750 case PS_TYPE_F64: \ 751 RESULT = VECTOR->data.F64[INDEX]; \ 752 break; \ 753 default: \ 754 RESULT = NAN; \ 755 break; \ 756 } 757 758 759 #define CHECK_LOWER_UPPER_BOUND(TABLE,INDEX,COLUMN) \ 760 switch (TABLE->index->type.type) { \ 761 case PS_TYPE_U8: \ 762 if( (psU8)index < TABLE->index->data.U8[0] ) { \ 763 *status = PS_LOOKUP_PAST_TOP; \ 764 } \ 765 if( (psU8)index > TABLE->index->data.U8[numRows-1]) { \ 766 *status = PS_LOOKUP_PAST_BOTTOM; \ 767 } \ 768 break; \ 769 case PS_TYPE_U16: \ 770 if( (psU16)index < TABLE->index->data.U16[0] ) { \ 771 *status = PS_LOOKUP_PAST_TOP; \ 772 } \ 773 if( (psU16)index > TABLE->index->data.U16[numRows-1] ) { \ 774 *status = PS_LOOKUP_PAST_BOTTOM; \ 775 } \ 776 break; \ 777 case PS_TYPE_U32: \ 778 if( (psU32)index < TABLE->index->data.U32[0]) { \ 779 *status = PS_LOOKUP_PAST_TOP; \ 780 } \ 781 if ( (psU32)index > TABLE->index->data.U32[numRows-1] ) { \ 782 *status = PS_LOOKUP_PAST_BOTTOM; \ 783 } \ 784 break; \ 785 case PS_TYPE_U64: \ 786 if( (psU64)index < TABLE->index->data.U64[0] ) { \ 787 *status = PS_LOOKUP_PAST_TOP; \ 788 } \ 789 if( (psU64)index > TABLE->index->data.U64[numRows-1] ) { \ 790 *status = PS_LOOKUP_PAST_BOTTOM; \ 791 } \ 792 break; \ 793 case PS_TYPE_S8: \ 794 if( (psS8)index < TABLE->index->data.S8[0] ) { \ 795 *status = PS_LOOKUP_PAST_TOP; \ 796 } \ 797 if( (psS8)index > TABLE->index->data.S8[numRows-1] ) { \ 798 *status = PS_LOOKUP_PAST_BOTTOM; \ 799 } \ 800 break; \ 801 case PS_TYPE_S16: \ 802 if( (psS16)index < TABLE->index->data.S16[0] ) { \ 803 *status = PS_LOOKUP_PAST_TOP; \ 804 } \ 805 if( (psS16)index > TABLE->index->data.S16[numRows-1] ) { \ 806 *status = PS_LOOKUP_PAST_BOTTOM; \ 807 } \ 808 break; \ 809 case PS_TYPE_S32: \ 810 if( (psS32)index < TABLE->index->data.S32[0] ) { \ 811 *status = PS_LOOKUP_PAST_TOP; \ 812 } \ 813 if( (psS32)index > TABLE->index->data.S32[numRows-1] ) { \ 814 *status = PS_LOOKUP_PAST_BOTTOM; \ 815 } \ 816 break; \ 817 case PS_TYPE_S64: \ 818 if( (psS64)index < TABLE->index->data.S64[0] ) { \ 819 *status = PS_LOOKUP_PAST_TOP; \ 820 } \ 821 if( (psS64)index > TABLE->index->data.S64[numRows-1] ) { \ 822 *status = PS_LOOKUP_PAST_BOTTOM; \ 823 } \ 824 break; \ 825 case PS_TYPE_F32: \ 826 if( (psF32)index < TABLE->index->data.F32[0] ) { \ 827 *status = PS_LOOKUP_PAST_TOP; \ 828 } \ 829 if( (psF32)index > TABLE->index->data.F32[numRows-1] ) { \ 830 *status = PS_LOOKUP_PAST_BOTTOM; \ 831 } \ 832 break; \ 833 case PS_TYPE_F64: \ 834 if( index < TABLE->index->data.F64[0] ) { \ 835 *status = PS_LOOKUP_PAST_TOP; \ 836 } \ 837 if( index > TABLE->index->data.F64[numRows-1] ) { \ 838 *status = PS_LOOKUP_PAST_BOTTOM; \ 839 } \ 840 break; \ 841 default: \ 842 *status = PS_LOOKUP_ERROR; \ 843 return NAN; \ 844 break; \ 845 } \ 846 if(*status == PS_LOOKUP_PAST_TOP) { \ 847 CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),0,out) \ 848 return out; \ 849 } else if (*status == PS_LOOKUP_PAST_BOTTOM) { \ 850 CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),numRows-1,out) \ 851 return out; \ 852 } 853 569 854 570 855 psF64 psLookupTableInterpolate(psLookupTable *table, psF64 index, psU64 column, psLookupStatusType *status) … … 576 861 psF64 out = 0.0; 577 862 psF64 denom = 0.0; 863 psF64 convertVal = 0.0; 864 psF64 tempVal = 0.0; 578 865 psVector *indexVec = NULL; 579 866 psVector *valuesVec = NULL; 580 867 psArray *values = NULL; 581 868 582 583 869 // Error checks 584 PS_PTR_CHECK_NULL(table,0.0); 870 // Set status to error prior to check since if checks fails it will immediately return 871 PS_PTR_CHECK_NULL(status,NAN); 872 *status = PS_LOOKUP_ERROR; 873 PS_PTR_CHECK_NULL(table,NAN); 585 874 indexVec = table->index; 586 875 values = table->values; 587 876 numRows = table->numRows; 588 877 numCols = table->numCols; 589 PS_PTR_CHECK_NULL(indexVec,0.0); 590 PS_PTR_CHECK_NULL(values,0.0); 591 PS_INT_CHECK_EQUALS(numRows, 0,0.0); 592 PS_INT_CHECK_EQUALS(numCols, 0,0.0); 593 PS_INT_CHECK_RANGE(column, 0, numCols-1, 0.0); 594 878 PS_PTR_CHECK_NULL(indexVec,NAN); 879 PS_PTR_CHECK_NULL(values,NAN); 880 PS_INT_CHECK_EQUALS(numRows, 0,NAN); 881 PS_INT_CHECK_EQUALS(numCols, 0,NAN); 882 PS_INT_CHECK_RANGE(column, 0, (int)(numCols-1), NAN); 595 883 596 884 valuesVec = (psVector*)values->data[column]; 597 PS_PTR_CHECK_NULL(indexVec,0.0); 598 599 if(index<indexVec->data.F64[0] || index<table->validFrom) { 600 601 // Index value past top of table 602 *status = PS_LOOKUP_PAST_TOP; 603 return valuesVec->data.F64[0]; 604 } else if(index>indexVec->data.F64[numRows-1] || index>table->validTo) { 605 606 // Index value past bottom of table 607 *status = PS_LOOKUP_PAST_BOTTOM; 608 return valuesVec->data.F64[numRows-1]; 885 PS_PTR_CHECK_NULL(indexVec,NAN); 886 887 // Set status to success since it passed all parameter checks 888 *status = PS_LOOKUP_SUCCESS; 889 890 // Verify the index is within the bounds of the table 891 CHECK_LOWER_UPPER_BOUND(table,index,column) 892 893 // Find location in table where specified index is between to entries 894 CONVERT_VALUE_TO_F64(indexVec, 0, convertVal) 895 while(index > convertVal ) { 896 hiIdx++; 897 if(hiIdx >= numRows) { 898 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 899 PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); 900 *status = PS_LOOKUP_ERROR; 901 return NAN; 902 } 903 CONVERT_VALUE_TO_F64(indexVec, hiIdx, convertVal) 904 } 905 906 // Check for negative low index and generate error 907 loIdx = hiIdx--; 908 if(loIdx < 0) { 909 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 910 PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); 911 *status = PS_LOOKUP_ERROR; 912 return NAN; 913 } 914 915 // Perform linear interpolation to calculate return value 916 CONVERT_VALUE_TO_F64(indexVec, hiIdx, denom) 917 CONVERT_VALUE_TO_F64(indexVec, loIdx, convertVal); 918 denom -= convertVal; 919 if(fabs(denom) < FLT_EPSILON) { 920 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); 921 *status = PS_LOOKUP_ERROR; 922 return NAN; 609 923 } else { 610 611 // Interpolation 612 while(index > indexVec->data.F64[hiIdx]) { 613 hiIdx++; 614 if(hiIdx >= numRows) { 615 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); 616 *status = PS_LOOKUP_ERROR; 617 return 0.0; 618 } 619 } 620 621 loIdx = hiIdx--; 622 if(loIdx < 0) { 623 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); 624 *status = PS_LOOKUP_ERROR; 625 return 0.0; 626 } 627 628 denom = indexVec->data.F64[hiIdx] - indexVec->data.F64[loIdx]; 629 if(fabs(denom) < FLT_EPSILON) { 630 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); 631 *status = PS_LOOKUP_ERROR; 632 return 0.0; 633 } else { 634 out = valuesVec->data.F64[loIdx]+(valuesVec->data.F64[hiIdx]-valuesVec->data.F64[loIdx]) 635 *(index-indexVec->data.F64[loIdx])/denom; 636 } 924 CONVERT_VALUE_TO_F64(valuesVec,hiIdx,tempVal) 925 CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) 926 tempVal -= convertVal; 927 CONVERT_VALUE_TO_F64(indexVec,loIdx,convertVal) 928 out = tempVal*(index-convertVal)/denom; 929 CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) 930 out += convertVal; 637 931 } 638 932 … … 646 940 psVector *outVector = NULL; 647 941 psLookupStatusType status = PS_LOOKUP_SUCCESS; 648 649 942 650 943 // Error checks … … 654 947 PS_INT_CHECK_EQUALS(numCols, 0,NULL); 655 948 656 outVector = psVectorAlloc(numCols , PS_TYPE_F64);949 outVector = psVectorAlloc(numCols+1, PS_TYPE_F64); 657 950 658 951 // Fill vectors with results and status of results -
trunk/psLib/src/dataIO/psLookupTable.h
r3264 r3271 7 7 * @author Ross Harman, MHPCC 8 8 * 9 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 * Copyright 2004 -2005Maui High Performance Computing Center, University of Hawaii12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 13 13 */ 14 14 -
trunk/psLib/src/fileUtils/psFileUtilsErrors.dat
r3028 r3271 15 15 psLookupTable_DIVIDE_BY_ZERO Divide by zero error during interpolation. 16 16 psLookupTable_INVALID_TYPE Invalid psLookupType, %d; 17 psLookupTable_TABLE_INVALID Lookup table is invalid. 17 18 # 18 19 psFits_NULL The input psFits object can not NULL. -
trunk/psLib/src/fileUtils/psFileUtilsErrors.h
r3264 r3271 7 7 * @author Robert DeSonia, MHPCC 8 8 * 9 * @version $Revision: 1.1 1$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 36 36 #define PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO "Divide by zero error during interpolation." 37 37 #define PS_ERRORTEXT_psLookupTable_INVALID_TYPE "Invalid psLookupType, %d;" 38 #define PS_ERRORTEXT_psLookupTable_TABLE_INVALID "Lookup table is invalid." 38 39 #define PS_ERRORTEXT_psFits_NULL "The input psFits object can not NULL." 39 40 #define PS_ERRORTEXT_psFits_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" -
trunk/psLib/src/fileUtils/psLookupTable.c
r3264 r3271 7 7 * @author Ross Harman, MHPCC 8 8 * 9 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 * Copyright 2004- 2005 Maui High Performance Computing Center, University of Hawaii12 * Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii 13 13 */ 14 14 #include <stdio.h> 15 15 #include <string.h> 16 16 #include <ctype.h> 17 //#ifdef DARWIN 17 18 #undef __STRICT_ANSI__ 19 //#endif 18 20 #include <stdlib.h> 19 #include <limits.h> 21 //#ifdef DARWIN 22 #define __STRICT_ANSI__ 23 //#endif 20 24 #include <math.h> 25 #include <stdlib.h> 21 26 22 27 #include "psMemory.h" … … 438 443 } 439 444 445 #define UPDATE_VALID_TO_FROM(TABLE) \ 446 switch (TABLE->index->type.type) { \ 447 case PS_TYPE_U8: \ 448 TABLE->validFrom = (psF64)TABLE->index->data.U8[0]; \ 449 TABLE->validTo = (psF64)TABLE->index->data.U8[TABLE->index->nalloc-1]; \ 450 break; \ 451 case PS_TYPE_S8: \ 452 TABLE->validFrom = (psF64)TABLE->index->data.S8[0]; \ 453 TABLE->validTo = (psF64)TABLE->index->data.S8[TABLE->index->nalloc-1]; \ 454 break; \ 455 case PS_TYPE_U16: \ 456 TABLE->validFrom = (psF64)TABLE->index->data.U16[0]; \ 457 TABLE->validTo = (psF64)TABLE->index->data.U16[TABLE->index->nalloc-1]; \ 458 break; \ 459 case PS_TYPE_S16: \ 460 TABLE->validFrom = (psF64)TABLE->index->data.S16[0]; \ 461 TABLE->validTo = (psF64)TABLE->index->data.S16[TABLE->index->nalloc-1]; \ 462 break; \ 463 case PS_TYPE_U32: \ 464 TABLE->validFrom = (psF64)TABLE->index->data.U32[0]; \ 465 TABLE->validTo = (psF64)TABLE->index->data.U32[TABLE->index->nalloc-1]; \ 466 break; \ 467 case PS_TYPE_S32: \ 468 TABLE->validFrom = (psF64)TABLE->index->data.S32[0]; \ 469 TABLE->validTo = (psF64)TABLE->index->data.S32[TABLE->index->nalloc-1]; \ 470 break; \ 471 case PS_TYPE_U64: \ 472 TABLE->validFrom = (psF64)TABLE->index->data.U64[0]; \ 473 TABLE->validTo = (psF64)TABLE->index->data.U64[TABLE->index->nalloc-1]; \ 474 break; \ 475 case PS_TYPE_S64: \ 476 TABLE->validFrom = (psF64)TABLE->index->data.S64[0]; \ 477 TABLE->validTo = (psF64)TABLE->index->data.S64[TABLE->index->nalloc-1]; \ 478 break; \ 479 case PS_TYPE_F32: \ 480 TABLE->validFrom = (psF64)TABLE->index->data.F32[0]; \ 481 TABLE->validTo = (psF64)TABLE->index->data.F32[TABLE->index->nalloc-1]; \ 482 break; \ 483 case PS_TYPE_F64: \ 484 TABLE->validFrom = (psF64)TABLE->index->data.F64[0]; \ 485 TABLE->validTo = (psF64)TABLE->index->data.F64[TABLE->index->nalloc-1]; \ 486 break; \ 487 default: \ 488 TABLE->validFrom = (psF64)0; \ 489 TABLE->validTo = (psF64)0; \ 490 break; \ 491 } 492 493 #define COPY_VECTOR_VALUES(VEC_OUT,INDEX_OUT,VEC_IN,INDEX_IN) \ 494 switch(((psVector*)(VEC_IN))->type.type) { \ 495 case PS_TYPE_U8: \ 496 ((psVector*)VEC_OUT)->data.U8[INDEX_OUT] = ((psVector*)VEC_IN)->data.U8[INDEX_IN]; \ 497 break; \ 498 case PS_TYPE_U16: \ 499 ((psVector*)VEC_OUT)->data.U16[INDEX_OUT] = ((psVector*)VEC_IN)->data.U16[INDEX_IN]; \ 500 break; \ 501 case PS_TYPE_U32: \ 502 ((psVector*)VEC_OUT)->data.U32[INDEX_OUT] = ((psVector*)VEC_IN)->data.U32[INDEX_IN]; \ 503 break; \ 504 case PS_TYPE_U64: \ 505 ((psVector*)VEC_OUT)->data.U64[INDEX_OUT] = ((psVector*)VEC_IN)->data.U64[INDEX_IN]; \ 506 break; \ 507 case PS_TYPE_S8: \ 508 ((psVector*)VEC_OUT)->data.S8[INDEX_OUT] = ((psVector*)VEC_IN)->data.S8[INDEX_IN]; \ 509 break; \ 510 case PS_TYPE_S16: \ 511 ((psVector*)VEC_OUT)->data.S16[INDEX_OUT] = ((psVector*)VEC_IN)->data.S16[INDEX_IN]; \ 512 break; \ 513 case PS_TYPE_S32: \ 514 ((psVector*)VEC_OUT)->data.S32[INDEX_OUT] = ((psVector*)VEC_IN)->data.S32[INDEX_IN]; \ 515 break; \ 516 case PS_TYPE_S64: \ 517 ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN]; \ 518 break; \ 519 case PS_TYPE_F32: \ 520 ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN]; \ 521 break; \ 522 case PS_TYPE_F64: \ 523 ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN]; \ 524 break; \ 525 default: \ 526 break; \ 527 } 528 529 440 530 psLookupTable* psLookupTableRead(psLookupTable *table) 441 531 { … … 446 536 char *linePtr = NULL; 447 537 psParseErrorType status = PS_PARSE_SUCCESS; 538 psParseErrorType lineStatus = PS_PARSE_SUCCESS; 448 539 psU64 lineCount = 0; 449 540 psU64 numRows = 0; … … 451 542 psU64 failedLines = 0; 452 543 FILE *fp = NULL; 453 psElemType elemType = PS_TYPE_PTR;544 psElemType elemType; 454 545 psVector *indexVec = NULL; 455 546 psVector *valuesVec = NULL; 456 547 psArray *values = NULL; 457 458 459 // Error checks548 psBool sortIndex = false; 549 550 // Check for NULL input table 460 551 PS_PTR_CHECK_NULL(table,NULL); 552 553 // Check for input table with NULL file name 461 554 PS_PTR_CHECK_NULL(table->fileName,NULL); 555 556 // Open table file specified by table->fileName 462 557 if((fp=fopen(table->fileName, "r")) == NULL) { 463 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, table->fileName); 558 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, 559 table->fileName); 464 560 return table; 465 561 } 466 562 563 // Initialize vector pointers 467 564 indexVec = table->index; 468 565 table->values = psArrayAlloc(1); … … 473 570 474 571 // Loop through file to get numRows, numCols, and column data types 475 while( fgets(line, MAX_STRING_LENGTH, fp) != NULL) {572 while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && (lineStatus == PS_PARSE_SUCCESS)) { 476 573 477 574 // Initialize variables for new line … … 485 582 486 583 // Determine column types from first line in data file after comments 487 while((strType=getToken(&linePtr, " ", &status)) != NULL) {584 while((strType=getToken(&linePtr," ",&status))) { 488 585 numCols++; 489 586 typeLine = false; … … 522 619 523 620 if(status) { 524 status =printError(lineCount, strValue, status);621 printError(lineCount, strValue, status); 525 622 failedLines++; 623 lineStatus = status; 624 status = PS_PARSE_SUCCESS; 526 625 } 527 626 psFree(strType); 528 627 } 529 628 } else { 530 531 629 // Parse and add values to all columns 532 630 numRows++; 533 631 numCols = 0; 534 while((strValue=getToken(&linePtr, " ", &status)) != NULL) {632 while((strValue=getToken(&linePtr," ",&status))) { 535 633 numCols++; 536 634 537 635 // Realloc number of rows as you go 538 636 if(numCols == 1) { 637 sortIndex = false; 539 638 indexVec = psVectorRecycle(indexVec, numRows, indexVec->type.type); 540 639 parseValue(indexVec, numRows-1, strValue, &status); … … 546 645 547 646 if(status) { 548 status =printError(lineCount, strValue, status);647 printError(lineCount, strValue, status); 549 648 failedLines++; 649 lineStatus = status; 650 status = PS_PARSE_SUCCESS; 550 651 } 551 652 psFree(strValue); … … 563 664 table->values = values; 564 665 666 // Set table values to indicate error detected during the read 667 if(lineStatus) { 668 table->numRows = 0; 669 table->numCols = 0; 670 table->validFrom = 0; 671 table->validTo = 0; 672 psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psLookupTable_TABLE_INVALID); 673 } else { 674 // Check if index vector needs to be sorted 675 psVector* sortedIndex = psVectorAlloc(table->numRows,PS_TYPE_U32); 676 sortedIndex = psVectorSortIndex(sortedIndex,table->index); 677 for(psS32 i = 0; i < numRows; i++ ) { 678 if(sortedIndex->data.U32[i] != i) { 679 sortIndex = true; 680 break; 681 } 682 } 683 // Check if it is necessary to sort value vectors 684 if(sortIndex) { 685 // Allocate new index vector 686 psVector* newIndexVector = psVectorAlloc(numRows,indexVec->type.type); 687 // Allocate new value vectors 688 psArray* newValueArray = psArrayAlloc(numCols-1); 689 for(psS32 j = 0; j < numCols-1; j++) { 690 psS32 type = ((psVector*)(table->values->data[j]))->type.type; 691 newValueArray->data[j] = psVectorAlloc(numRows,type); 692 } 693 for(psS32 i = 0; i < numRows; i++) { 694 // Populate new index vector 695 psU32 sortIndex = sortedIndex->data.U32[i]; 696 COPY_VECTOR_VALUES(newIndexVector,i, indexVec,sortIndex) 697 // For every column populate new value vectors 698 for(psS32 j=0; j < numCols-1; j++) { 699 COPY_VECTOR_VALUES(newValueArray->data[j],i,table->values->data[j], sortIndex) 700 } 701 } 702 // Free old index vector 703 psFree(table->index); 704 // Free old value vectors 705 psFree(table->values); 706 // Assign new vector value to table 707 table->index = newIndexVector; 708 // Assign new value vectors to table array 709 table->values = newValueArray; 710 } 711 psFree(sortedIndex); 712 // Update the validTo and validFrom 713 UPDATE_VALID_TO_FROM(table) 714 } 715 565 716 fclose(fp); 566 717 567 718 return table; 568 719 } 720 721 #define CONVERT_VALUE_TO_F64(VECTOR,INDEX,RESULT) \ 722 switch(VECTOR->type.type) { \ 723 case PS_TYPE_U8: \ 724 RESULT = (psF64)VECTOR->data.U8[INDEX]; \ 725 break; \ 726 case PS_TYPE_U16: \ 727 RESULT = (psF64)VECTOR->data.U16[INDEX]; \ 728 break; \ 729 case PS_TYPE_U32: \ 730 RESULT = (psF64)VECTOR->data.U32[INDEX]; \ 731 break; \ 732 case PS_TYPE_U64: \ 733 RESULT = (psF64)VECTOR->data.U64[INDEX]; \ 734 break; \ 735 case PS_TYPE_S8: \ 736 RESULT = (psF64)VECTOR->data.S8[INDEX]; \ 737 break; \ 738 case PS_TYPE_S16: \ 739 RESULT = (psF64)VECTOR->data.S16[INDEX]; \ 740 break; \ 741 case PS_TYPE_S32: \ 742 RESULT = (psF64)VECTOR->data.S32[INDEX]; \ 743 break; \ 744 case PS_TYPE_S64: \ 745 RESULT = (psF64)VECTOR->data.S64[INDEX]; \ 746 break; \ 747 case PS_TYPE_F32: \ 748 RESULT = (psF64)VECTOR->data.F32[INDEX]; \ 749 break; \ 750 case PS_TYPE_F64: \ 751 RESULT = VECTOR->data.F64[INDEX]; \ 752 break; \ 753 default: \ 754 RESULT = NAN; \ 755 break; \ 756 } 757 758 759 #define CHECK_LOWER_UPPER_BOUND(TABLE,INDEX,COLUMN) \ 760 switch (TABLE->index->type.type) { \ 761 case PS_TYPE_U8: \ 762 if( (psU8)index < TABLE->index->data.U8[0] ) { \ 763 *status = PS_LOOKUP_PAST_TOP; \ 764 } \ 765 if( (psU8)index > TABLE->index->data.U8[numRows-1]) { \ 766 *status = PS_LOOKUP_PAST_BOTTOM; \ 767 } \ 768 break; \ 769 case PS_TYPE_U16: \ 770 if( (psU16)index < TABLE->index->data.U16[0] ) { \ 771 *status = PS_LOOKUP_PAST_TOP; \ 772 } \ 773 if( (psU16)index > TABLE->index->data.U16[numRows-1] ) { \ 774 *status = PS_LOOKUP_PAST_BOTTOM; \ 775 } \ 776 break; \ 777 case PS_TYPE_U32: \ 778 if( (psU32)index < TABLE->index->data.U32[0]) { \ 779 *status = PS_LOOKUP_PAST_TOP; \ 780 } \ 781 if ( (psU32)index > TABLE->index->data.U32[numRows-1] ) { \ 782 *status = PS_LOOKUP_PAST_BOTTOM; \ 783 } \ 784 break; \ 785 case PS_TYPE_U64: \ 786 if( (psU64)index < TABLE->index->data.U64[0] ) { \ 787 *status = PS_LOOKUP_PAST_TOP; \ 788 } \ 789 if( (psU64)index > TABLE->index->data.U64[numRows-1] ) { \ 790 *status = PS_LOOKUP_PAST_BOTTOM; \ 791 } \ 792 break; \ 793 case PS_TYPE_S8: \ 794 if( (psS8)index < TABLE->index->data.S8[0] ) { \ 795 *status = PS_LOOKUP_PAST_TOP; \ 796 } \ 797 if( (psS8)index > TABLE->index->data.S8[numRows-1] ) { \ 798 *status = PS_LOOKUP_PAST_BOTTOM; \ 799 } \ 800 break; \ 801 case PS_TYPE_S16: \ 802 if( (psS16)index < TABLE->index->data.S16[0] ) { \ 803 *status = PS_LOOKUP_PAST_TOP; \ 804 } \ 805 if( (psS16)index > TABLE->index->data.S16[numRows-1] ) { \ 806 *status = PS_LOOKUP_PAST_BOTTOM; \ 807 } \ 808 break; \ 809 case PS_TYPE_S32: \ 810 if( (psS32)index < TABLE->index->data.S32[0] ) { \ 811 *status = PS_LOOKUP_PAST_TOP; \ 812 } \ 813 if( (psS32)index > TABLE->index->data.S32[numRows-1] ) { \ 814 *status = PS_LOOKUP_PAST_BOTTOM; \ 815 } \ 816 break; \ 817 case PS_TYPE_S64: \ 818 if( (psS64)index < TABLE->index->data.S64[0] ) { \ 819 *status = PS_LOOKUP_PAST_TOP; \ 820 } \ 821 if( (psS64)index > TABLE->index->data.S64[numRows-1] ) { \ 822 *status = PS_LOOKUP_PAST_BOTTOM; \ 823 } \ 824 break; \ 825 case PS_TYPE_F32: \ 826 if( (psF32)index < TABLE->index->data.F32[0] ) { \ 827 *status = PS_LOOKUP_PAST_TOP; \ 828 } \ 829 if( (psF32)index > TABLE->index->data.F32[numRows-1] ) { \ 830 *status = PS_LOOKUP_PAST_BOTTOM; \ 831 } \ 832 break; \ 833 case PS_TYPE_F64: \ 834 if( index < TABLE->index->data.F64[0] ) { \ 835 *status = PS_LOOKUP_PAST_TOP; \ 836 } \ 837 if( index > TABLE->index->data.F64[numRows-1] ) { \ 838 *status = PS_LOOKUP_PAST_BOTTOM; \ 839 } \ 840 break; \ 841 default: \ 842 *status = PS_LOOKUP_ERROR; \ 843 return NAN; \ 844 break; \ 845 } \ 846 if(*status == PS_LOOKUP_PAST_TOP) { \ 847 CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),0,out) \ 848 return out; \ 849 } else if (*status == PS_LOOKUP_PAST_BOTTOM) { \ 850 CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),numRows-1,out) \ 851 return out; \ 852 } 853 569 854 570 855 psF64 psLookupTableInterpolate(psLookupTable *table, psF64 index, psU64 column, psLookupStatusType *status) … … 576 861 psF64 out = 0.0; 577 862 psF64 denom = 0.0; 863 psF64 convertVal = 0.0; 864 psF64 tempVal = 0.0; 578 865 psVector *indexVec = NULL; 579 866 psVector *valuesVec = NULL; 580 867 psArray *values = NULL; 581 868 582 583 869 // Error checks 584 PS_PTR_CHECK_NULL(table,0.0); 870 // Set status to error prior to check since if checks fails it will immediately return 871 PS_PTR_CHECK_NULL(status,NAN); 872 *status = PS_LOOKUP_ERROR; 873 PS_PTR_CHECK_NULL(table,NAN); 585 874 indexVec = table->index; 586 875 values = table->values; 587 876 numRows = table->numRows; 588 877 numCols = table->numCols; 589 PS_PTR_CHECK_NULL(indexVec,0.0); 590 PS_PTR_CHECK_NULL(values,0.0); 591 PS_INT_CHECK_EQUALS(numRows, 0,0.0); 592 PS_INT_CHECK_EQUALS(numCols, 0,0.0); 593 PS_INT_CHECK_RANGE(column, 0, numCols-1, 0.0); 594 878 PS_PTR_CHECK_NULL(indexVec,NAN); 879 PS_PTR_CHECK_NULL(values,NAN); 880 PS_INT_CHECK_EQUALS(numRows, 0,NAN); 881 PS_INT_CHECK_EQUALS(numCols, 0,NAN); 882 PS_INT_CHECK_RANGE(column, 0, (int)(numCols-1), NAN); 595 883 596 884 valuesVec = (psVector*)values->data[column]; 597 PS_PTR_CHECK_NULL(indexVec,0.0); 598 599 if(index<indexVec->data.F64[0] || index<table->validFrom) { 600 601 // Index value past top of table 602 *status = PS_LOOKUP_PAST_TOP; 603 return valuesVec->data.F64[0]; 604 } else if(index>indexVec->data.F64[numRows-1] || index>table->validTo) { 605 606 // Index value past bottom of table 607 *status = PS_LOOKUP_PAST_BOTTOM; 608 return valuesVec->data.F64[numRows-1]; 885 PS_PTR_CHECK_NULL(indexVec,NAN); 886 887 // Set status to success since it passed all parameter checks 888 *status = PS_LOOKUP_SUCCESS; 889 890 // Verify the index is within the bounds of the table 891 CHECK_LOWER_UPPER_BOUND(table,index,column) 892 893 // Find location in table where specified index is between to entries 894 CONVERT_VALUE_TO_F64(indexVec, 0, convertVal) 895 while(index > convertVal ) { 896 hiIdx++; 897 if(hiIdx >= numRows) { 898 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 899 PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); 900 *status = PS_LOOKUP_ERROR; 901 return NAN; 902 } 903 CONVERT_VALUE_TO_F64(indexVec, hiIdx, convertVal) 904 } 905 906 // Check for negative low index and generate error 907 loIdx = hiIdx--; 908 if(loIdx < 0) { 909 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 910 PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); 911 *status = PS_LOOKUP_ERROR; 912 return NAN; 913 } 914 915 // Perform linear interpolation to calculate return value 916 CONVERT_VALUE_TO_F64(indexVec, hiIdx, denom) 917 CONVERT_VALUE_TO_F64(indexVec, loIdx, convertVal); 918 denom -= convertVal; 919 if(fabs(denom) < FLT_EPSILON) { 920 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); 921 *status = PS_LOOKUP_ERROR; 922 return NAN; 609 923 } else { 610 611 // Interpolation 612 while(index > indexVec->data.F64[hiIdx]) { 613 hiIdx++; 614 if(hiIdx >= numRows) { 615 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); 616 *status = PS_LOOKUP_ERROR; 617 return 0.0; 618 } 619 } 620 621 loIdx = hiIdx--; 622 if(loIdx < 0) { 623 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); 624 *status = PS_LOOKUP_ERROR; 625 return 0.0; 626 } 627 628 denom = indexVec->data.F64[hiIdx] - indexVec->data.F64[loIdx]; 629 if(fabs(denom) < FLT_EPSILON) { 630 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); 631 *status = PS_LOOKUP_ERROR; 632 return 0.0; 633 } else { 634 out = valuesVec->data.F64[loIdx]+(valuesVec->data.F64[hiIdx]-valuesVec->data.F64[loIdx]) 635 *(index-indexVec->data.F64[loIdx])/denom; 636 } 924 CONVERT_VALUE_TO_F64(valuesVec,hiIdx,tempVal) 925 CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) 926 tempVal -= convertVal; 927 CONVERT_VALUE_TO_F64(indexVec,loIdx,convertVal) 928 out = tempVal*(index-convertVal)/denom; 929 CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) 930 out += convertVal; 637 931 } 638 932 … … 646 940 psVector *outVector = NULL; 647 941 psLookupStatusType status = PS_LOOKUP_SUCCESS; 648 649 942 650 943 // Error checks … … 654 947 PS_INT_CHECK_EQUALS(numCols, 0,NULL); 655 948 656 outVector = psVectorAlloc(numCols , PS_TYPE_F64);949 outVector = psVectorAlloc(numCols+1, PS_TYPE_F64); 657 950 658 951 // Fill vectors with results and status of results -
trunk/psLib/src/fileUtils/psLookupTable.h
r3264 r3271 7 7 * @author Ross Harman, MHPCC 8 8 * 9 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 * Copyright 2004 -2005Maui High Performance Computing Center, University of Hawaii12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 13 13 */ 14 14 -
trunk/psLib/src/types/psLookupTable.c
r3264 r3271 7 7 * @author Ross Harman, MHPCC 8 8 * 9 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 * Copyright 2004- 2005 Maui High Performance Computing Center, University of Hawaii12 * Copyright 2004-5 Maui High Performance Computing Center, University of Hawaii 13 13 */ 14 14 #include <stdio.h> 15 15 #include <string.h> 16 16 #include <ctype.h> 17 //#ifdef DARWIN 17 18 #undef __STRICT_ANSI__ 19 //#endif 18 20 #include <stdlib.h> 19 #include <limits.h> 21 //#ifdef DARWIN 22 #define __STRICT_ANSI__ 23 //#endif 20 24 #include <math.h> 25 #include <stdlib.h> 21 26 22 27 #include "psMemory.h" … … 438 443 } 439 444 445 #define UPDATE_VALID_TO_FROM(TABLE) \ 446 switch (TABLE->index->type.type) { \ 447 case PS_TYPE_U8: \ 448 TABLE->validFrom = (psF64)TABLE->index->data.U8[0]; \ 449 TABLE->validTo = (psF64)TABLE->index->data.U8[TABLE->index->nalloc-1]; \ 450 break; \ 451 case PS_TYPE_S8: \ 452 TABLE->validFrom = (psF64)TABLE->index->data.S8[0]; \ 453 TABLE->validTo = (psF64)TABLE->index->data.S8[TABLE->index->nalloc-1]; \ 454 break; \ 455 case PS_TYPE_U16: \ 456 TABLE->validFrom = (psF64)TABLE->index->data.U16[0]; \ 457 TABLE->validTo = (psF64)TABLE->index->data.U16[TABLE->index->nalloc-1]; \ 458 break; \ 459 case PS_TYPE_S16: \ 460 TABLE->validFrom = (psF64)TABLE->index->data.S16[0]; \ 461 TABLE->validTo = (psF64)TABLE->index->data.S16[TABLE->index->nalloc-1]; \ 462 break; \ 463 case PS_TYPE_U32: \ 464 TABLE->validFrom = (psF64)TABLE->index->data.U32[0]; \ 465 TABLE->validTo = (psF64)TABLE->index->data.U32[TABLE->index->nalloc-1]; \ 466 break; \ 467 case PS_TYPE_S32: \ 468 TABLE->validFrom = (psF64)TABLE->index->data.S32[0]; \ 469 TABLE->validTo = (psF64)TABLE->index->data.S32[TABLE->index->nalloc-1]; \ 470 break; \ 471 case PS_TYPE_U64: \ 472 TABLE->validFrom = (psF64)TABLE->index->data.U64[0]; \ 473 TABLE->validTo = (psF64)TABLE->index->data.U64[TABLE->index->nalloc-1]; \ 474 break; \ 475 case PS_TYPE_S64: \ 476 TABLE->validFrom = (psF64)TABLE->index->data.S64[0]; \ 477 TABLE->validTo = (psF64)TABLE->index->data.S64[TABLE->index->nalloc-1]; \ 478 break; \ 479 case PS_TYPE_F32: \ 480 TABLE->validFrom = (psF64)TABLE->index->data.F32[0]; \ 481 TABLE->validTo = (psF64)TABLE->index->data.F32[TABLE->index->nalloc-1]; \ 482 break; \ 483 case PS_TYPE_F64: \ 484 TABLE->validFrom = (psF64)TABLE->index->data.F64[0]; \ 485 TABLE->validTo = (psF64)TABLE->index->data.F64[TABLE->index->nalloc-1]; \ 486 break; \ 487 default: \ 488 TABLE->validFrom = (psF64)0; \ 489 TABLE->validTo = (psF64)0; \ 490 break; \ 491 } 492 493 #define COPY_VECTOR_VALUES(VEC_OUT,INDEX_OUT,VEC_IN,INDEX_IN) \ 494 switch(((psVector*)(VEC_IN))->type.type) { \ 495 case PS_TYPE_U8: \ 496 ((psVector*)VEC_OUT)->data.U8[INDEX_OUT] = ((psVector*)VEC_IN)->data.U8[INDEX_IN]; \ 497 break; \ 498 case PS_TYPE_U16: \ 499 ((psVector*)VEC_OUT)->data.U16[INDEX_OUT] = ((psVector*)VEC_IN)->data.U16[INDEX_IN]; \ 500 break; \ 501 case PS_TYPE_U32: \ 502 ((psVector*)VEC_OUT)->data.U32[INDEX_OUT] = ((psVector*)VEC_IN)->data.U32[INDEX_IN]; \ 503 break; \ 504 case PS_TYPE_U64: \ 505 ((psVector*)VEC_OUT)->data.U64[INDEX_OUT] = ((psVector*)VEC_IN)->data.U64[INDEX_IN]; \ 506 break; \ 507 case PS_TYPE_S8: \ 508 ((psVector*)VEC_OUT)->data.S8[INDEX_OUT] = ((psVector*)VEC_IN)->data.S8[INDEX_IN]; \ 509 break; \ 510 case PS_TYPE_S16: \ 511 ((psVector*)VEC_OUT)->data.S16[INDEX_OUT] = ((psVector*)VEC_IN)->data.S16[INDEX_IN]; \ 512 break; \ 513 case PS_TYPE_S32: \ 514 ((psVector*)VEC_OUT)->data.S32[INDEX_OUT] = ((psVector*)VEC_IN)->data.S32[INDEX_IN]; \ 515 break; \ 516 case PS_TYPE_S64: \ 517 ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN]; \ 518 break; \ 519 case PS_TYPE_F32: \ 520 ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN]; \ 521 break; \ 522 case PS_TYPE_F64: \ 523 ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN]; \ 524 break; \ 525 default: \ 526 break; \ 527 } 528 529 440 530 psLookupTable* psLookupTableRead(psLookupTable *table) 441 531 { … … 446 536 char *linePtr = NULL; 447 537 psParseErrorType status = PS_PARSE_SUCCESS; 538 psParseErrorType lineStatus = PS_PARSE_SUCCESS; 448 539 psU64 lineCount = 0; 449 540 psU64 numRows = 0; … … 451 542 psU64 failedLines = 0; 452 543 FILE *fp = NULL; 453 psElemType elemType = PS_TYPE_PTR;544 psElemType elemType; 454 545 psVector *indexVec = NULL; 455 546 psVector *valuesVec = NULL; 456 547 psArray *values = NULL; 457 458 459 // Error checks548 psBool sortIndex = false; 549 550 // Check for NULL input table 460 551 PS_PTR_CHECK_NULL(table,NULL); 552 553 // Check for input table with NULL file name 461 554 PS_PTR_CHECK_NULL(table->fileName,NULL); 555 556 // Open table file specified by table->fileName 462 557 if((fp=fopen(table->fileName, "r")) == NULL) { 463 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, table->fileName); 558 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_FILE_NOT_FOUND, 559 table->fileName); 464 560 return table; 465 561 } 466 562 563 // Initialize vector pointers 467 564 indexVec = table->index; 468 565 table->values = psArrayAlloc(1); … … 473 570 474 571 // Loop through file to get numRows, numCols, and column data types 475 while( fgets(line, MAX_STRING_LENGTH, fp) != NULL) {572 while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) && (lineStatus == PS_PARSE_SUCCESS)) { 476 573 477 574 // Initialize variables for new line … … 485 582 486 583 // Determine column types from first line in data file after comments 487 while((strType=getToken(&linePtr, " ", &status)) != NULL) {584 while((strType=getToken(&linePtr," ",&status))) { 488 585 numCols++; 489 586 typeLine = false; … … 522 619 523 620 if(status) { 524 status =printError(lineCount, strValue, status);621 printError(lineCount, strValue, status); 525 622 failedLines++; 623 lineStatus = status; 624 status = PS_PARSE_SUCCESS; 526 625 } 527 626 psFree(strType); 528 627 } 529 628 } else { 530 531 629 // Parse and add values to all columns 532 630 numRows++; 533 631 numCols = 0; 534 while((strValue=getToken(&linePtr, " ", &status)) != NULL) {632 while((strValue=getToken(&linePtr," ",&status))) { 535 633 numCols++; 536 634 537 635 // Realloc number of rows as you go 538 636 if(numCols == 1) { 637 sortIndex = false; 539 638 indexVec = psVectorRecycle(indexVec, numRows, indexVec->type.type); 540 639 parseValue(indexVec, numRows-1, strValue, &status); … … 546 645 547 646 if(status) { 548 status =printError(lineCount, strValue, status);647 printError(lineCount, strValue, status); 549 648 failedLines++; 649 lineStatus = status; 650 status = PS_PARSE_SUCCESS; 550 651 } 551 652 psFree(strValue); … … 563 664 table->values = values; 564 665 666 // Set table values to indicate error detected during the read 667 if(lineStatus) { 668 table->numRows = 0; 669 table->numCols = 0; 670 table->validFrom = 0; 671 table->validTo = 0; 672 psError(PS_ERR_UNKNOWN,false,PS_ERRORTEXT_psLookupTable_TABLE_INVALID); 673 } else { 674 // Check if index vector needs to be sorted 675 psVector* sortedIndex = psVectorAlloc(table->numRows,PS_TYPE_U32); 676 sortedIndex = psVectorSortIndex(sortedIndex,table->index); 677 for(psS32 i = 0; i < numRows; i++ ) { 678 if(sortedIndex->data.U32[i] != i) { 679 sortIndex = true; 680 break; 681 } 682 } 683 // Check if it is necessary to sort value vectors 684 if(sortIndex) { 685 // Allocate new index vector 686 psVector* newIndexVector = psVectorAlloc(numRows,indexVec->type.type); 687 // Allocate new value vectors 688 psArray* newValueArray = psArrayAlloc(numCols-1); 689 for(psS32 j = 0; j < numCols-1; j++) { 690 psS32 type = ((psVector*)(table->values->data[j]))->type.type; 691 newValueArray->data[j] = psVectorAlloc(numRows,type); 692 } 693 for(psS32 i = 0; i < numRows; i++) { 694 // Populate new index vector 695 psU32 sortIndex = sortedIndex->data.U32[i]; 696 COPY_VECTOR_VALUES(newIndexVector,i, indexVec,sortIndex) 697 // For every column populate new value vectors 698 for(psS32 j=0; j < numCols-1; j++) { 699 COPY_VECTOR_VALUES(newValueArray->data[j],i,table->values->data[j], sortIndex) 700 } 701 } 702 // Free old index vector 703 psFree(table->index); 704 // Free old value vectors 705 psFree(table->values); 706 // Assign new vector value to table 707 table->index = newIndexVector; 708 // Assign new value vectors to table array 709 table->values = newValueArray; 710 } 711 psFree(sortedIndex); 712 // Update the validTo and validFrom 713 UPDATE_VALID_TO_FROM(table) 714 } 715 565 716 fclose(fp); 566 717 567 718 return table; 568 719 } 720 721 #define CONVERT_VALUE_TO_F64(VECTOR,INDEX,RESULT) \ 722 switch(VECTOR->type.type) { \ 723 case PS_TYPE_U8: \ 724 RESULT = (psF64)VECTOR->data.U8[INDEX]; \ 725 break; \ 726 case PS_TYPE_U16: \ 727 RESULT = (psF64)VECTOR->data.U16[INDEX]; \ 728 break; \ 729 case PS_TYPE_U32: \ 730 RESULT = (psF64)VECTOR->data.U32[INDEX]; \ 731 break; \ 732 case PS_TYPE_U64: \ 733 RESULT = (psF64)VECTOR->data.U64[INDEX]; \ 734 break; \ 735 case PS_TYPE_S8: \ 736 RESULT = (psF64)VECTOR->data.S8[INDEX]; \ 737 break; \ 738 case PS_TYPE_S16: \ 739 RESULT = (psF64)VECTOR->data.S16[INDEX]; \ 740 break; \ 741 case PS_TYPE_S32: \ 742 RESULT = (psF64)VECTOR->data.S32[INDEX]; \ 743 break; \ 744 case PS_TYPE_S64: \ 745 RESULT = (psF64)VECTOR->data.S64[INDEX]; \ 746 break; \ 747 case PS_TYPE_F32: \ 748 RESULT = (psF64)VECTOR->data.F32[INDEX]; \ 749 break; \ 750 case PS_TYPE_F64: \ 751 RESULT = VECTOR->data.F64[INDEX]; \ 752 break; \ 753 default: \ 754 RESULT = NAN; \ 755 break; \ 756 } 757 758 759 #define CHECK_LOWER_UPPER_BOUND(TABLE,INDEX,COLUMN) \ 760 switch (TABLE->index->type.type) { \ 761 case PS_TYPE_U8: \ 762 if( (psU8)index < TABLE->index->data.U8[0] ) { \ 763 *status = PS_LOOKUP_PAST_TOP; \ 764 } \ 765 if( (psU8)index > TABLE->index->data.U8[numRows-1]) { \ 766 *status = PS_LOOKUP_PAST_BOTTOM; \ 767 } \ 768 break; \ 769 case PS_TYPE_U16: \ 770 if( (psU16)index < TABLE->index->data.U16[0] ) { \ 771 *status = PS_LOOKUP_PAST_TOP; \ 772 } \ 773 if( (psU16)index > TABLE->index->data.U16[numRows-1] ) { \ 774 *status = PS_LOOKUP_PAST_BOTTOM; \ 775 } \ 776 break; \ 777 case PS_TYPE_U32: \ 778 if( (psU32)index < TABLE->index->data.U32[0]) { \ 779 *status = PS_LOOKUP_PAST_TOP; \ 780 } \ 781 if ( (psU32)index > TABLE->index->data.U32[numRows-1] ) { \ 782 *status = PS_LOOKUP_PAST_BOTTOM; \ 783 } \ 784 break; \ 785 case PS_TYPE_U64: \ 786 if( (psU64)index < TABLE->index->data.U64[0] ) { \ 787 *status = PS_LOOKUP_PAST_TOP; \ 788 } \ 789 if( (psU64)index > TABLE->index->data.U64[numRows-1] ) { \ 790 *status = PS_LOOKUP_PAST_BOTTOM; \ 791 } \ 792 break; \ 793 case PS_TYPE_S8: \ 794 if( (psS8)index < TABLE->index->data.S8[0] ) { \ 795 *status = PS_LOOKUP_PAST_TOP; \ 796 } \ 797 if( (psS8)index > TABLE->index->data.S8[numRows-1] ) { \ 798 *status = PS_LOOKUP_PAST_BOTTOM; \ 799 } \ 800 break; \ 801 case PS_TYPE_S16: \ 802 if( (psS16)index < TABLE->index->data.S16[0] ) { \ 803 *status = PS_LOOKUP_PAST_TOP; \ 804 } \ 805 if( (psS16)index > TABLE->index->data.S16[numRows-1] ) { \ 806 *status = PS_LOOKUP_PAST_BOTTOM; \ 807 } \ 808 break; \ 809 case PS_TYPE_S32: \ 810 if( (psS32)index < TABLE->index->data.S32[0] ) { \ 811 *status = PS_LOOKUP_PAST_TOP; \ 812 } \ 813 if( (psS32)index > TABLE->index->data.S32[numRows-1] ) { \ 814 *status = PS_LOOKUP_PAST_BOTTOM; \ 815 } \ 816 break; \ 817 case PS_TYPE_S64: \ 818 if( (psS64)index < TABLE->index->data.S64[0] ) { \ 819 *status = PS_LOOKUP_PAST_TOP; \ 820 } \ 821 if( (psS64)index > TABLE->index->data.S64[numRows-1] ) { \ 822 *status = PS_LOOKUP_PAST_BOTTOM; \ 823 } \ 824 break; \ 825 case PS_TYPE_F32: \ 826 if( (psF32)index < TABLE->index->data.F32[0] ) { \ 827 *status = PS_LOOKUP_PAST_TOP; \ 828 } \ 829 if( (psF32)index > TABLE->index->data.F32[numRows-1] ) { \ 830 *status = PS_LOOKUP_PAST_BOTTOM; \ 831 } \ 832 break; \ 833 case PS_TYPE_F64: \ 834 if( index < TABLE->index->data.F64[0] ) { \ 835 *status = PS_LOOKUP_PAST_TOP; \ 836 } \ 837 if( index > TABLE->index->data.F64[numRows-1] ) { \ 838 *status = PS_LOOKUP_PAST_BOTTOM; \ 839 } \ 840 break; \ 841 default: \ 842 *status = PS_LOOKUP_ERROR; \ 843 return NAN; \ 844 break; \ 845 } \ 846 if(*status == PS_LOOKUP_PAST_TOP) { \ 847 CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),0,out) \ 848 return out; \ 849 } else if (*status == PS_LOOKUP_PAST_BOTTOM) { \ 850 CONVERT_VALUE_TO_F64(((psVector*)(TABLE->values->data[COLUMN])),numRows-1,out) \ 851 return out; \ 852 } 853 569 854 570 855 psF64 psLookupTableInterpolate(psLookupTable *table, psF64 index, psU64 column, psLookupStatusType *status) … … 576 861 psF64 out = 0.0; 577 862 psF64 denom = 0.0; 863 psF64 convertVal = 0.0; 864 psF64 tempVal = 0.0; 578 865 psVector *indexVec = NULL; 579 866 psVector *valuesVec = NULL; 580 867 psArray *values = NULL; 581 868 582 583 869 // Error checks 584 PS_PTR_CHECK_NULL(table,0.0); 870 // Set status to error prior to check since if checks fails it will immediately return 871 PS_PTR_CHECK_NULL(status,NAN); 872 *status = PS_LOOKUP_ERROR; 873 PS_PTR_CHECK_NULL(table,NAN); 585 874 indexVec = table->index; 586 875 values = table->values; 587 876 numRows = table->numRows; 588 877 numCols = table->numCols; 589 PS_PTR_CHECK_NULL(indexVec,0.0); 590 PS_PTR_CHECK_NULL(values,0.0); 591 PS_INT_CHECK_EQUALS(numRows, 0,0.0); 592 PS_INT_CHECK_EQUALS(numCols, 0,0.0); 593 PS_INT_CHECK_RANGE(column, 0, numCols-1, 0.0); 594 878 PS_PTR_CHECK_NULL(indexVec,NAN); 879 PS_PTR_CHECK_NULL(values,NAN); 880 PS_INT_CHECK_EQUALS(numRows, 0,NAN); 881 PS_INT_CHECK_EQUALS(numCols, 0,NAN); 882 PS_INT_CHECK_RANGE(column, 0, (int)(numCols-1), NAN); 595 883 596 884 valuesVec = (psVector*)values->data[column]; 597 PS_PTR_CHECK_NULL(indexVec,0.0); 598 599 if(index<indexVec->data.F64[0] || index<table->validFrom) { 600 601 // Index value past top of table 602 *status = PS_LOOKUP_PAST_TOP; 603 return valuesVec->data.F64[0]; 604 } else if(index>indexVec->data.F64[numRows-1] || index>table->validTo) { 605 606 // Index value past bottom of table 607 *status = PS_LOOKUP_PAST_BOTTOM; 608 return valuesVec->data.F64[numRows-1]; 885 PS_PTR_CHECK_NULL(indexVec,NAN); 886 887 // Set status to success since it passed all parameter checks 888 *status = PS_LOOKUP_SUCCESS; 889 890 // Verify the index is within the bounds of the table 891 CHECK_LOWER_UPPER_BOUND(table,index,column) 892 893 // Find location in table where specified index is between to entries 894 CONVERT_VALUE_TO_F64(indexVec, 0, convertVal) 895 while(index > convertVal ) { 896 hiIdx++; 897 if(hiIdx >= numRows) { 898 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 899 PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); 900 *status = PS_LOOKUP_ERROR; 901 return NAN; 902 } 903 CONVERT_VALUE_TO_F64(indexVec, hiIdx, convertVal) 904 } 905 906 // Check for negative low index and generate error 907 loIdx = hiIdx--; 908 if(loIdx < 0) { 909 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 910 PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); 911 *status = PS_LOOKUP_ERROR; 912 return NAN; 913 } 914 915 // Perform linear interpolation to calculate return value 916 CONVERT_VALUE_TO_F64(indexVec, hiIdx, denom) 917 CONVERT_VALUE_TO_F64(indexVec, loIdx, convertVal); 918 denom -= convertVal; 919 if(fabs(denom) < FLT_EPSILON) { 920 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); 921 *status = PS_LOOKUP_ERROR; 922 return NAN; 609 923 } else { 610 611 // Interpolation 612 while(index > indexVec->data.F64[hiIdx]) { 613 hiIdx++; 614 if(hiIdx >= numRows) { 615 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_HIGH, hiIdx); 616 *status = PS_LOOKUP_ERROR; 617 return 0.0; 618 } 619 } 620 621 loIdx = hiIdx--; 622 if(loIdx < 0) { 623 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_INTERPOLATE_LOW, loIdx); 624 *status = PS_LOOKUP_ERROR; 625 return 0.0; 626 } 627 628 denom = indexVec->data.F64[hiIdx] - indexVec->data.F64[loIdx]; 629 if(fabs(denom) < FLT_EPSILON) { 630 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psLookupTable_DIVIDE_BY_ZERO); 631 *status = PS_LOOKUP_ERROR; 632 return 0.0; 633 } else { 634 out = valuesVec->data.F64[loIdx]+(valuesVec->data.F64[hiIdx]-valuesVec->data.F64[loIdx]) 635 *(index-indexVec->data.F64[loIdx])/denom; 636 } 924 CONVERT_VALUE_TO_F64(valuesVec,hiIdx,tempVal) 925 CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) 926 tempVal -= convertVal; 927 CONVERT_VALUE_TO_F64(indexVec,loIdx,convertVal) 928 out = tempVal*(index-convertVal)/denom; 929 CONVERT_VALUE_TO_F64(valuesVec,loIdx,convertVal) 930 out += convertVal; 637 931 } 638 932 … … 646 940 psVector *outVector = NULL; 647 941 psLookupStatusType status = PS_LOOKUP_SUCCESS; 648 649 942 650 943 // Error checks … … 654 947 PS_INT_CHECK_EQUALS(numCols, 0,NULL); 655 948 656 outVector = psVectorAlloc(numCols , PS_TYPE_F64);949 outVector = psVectorAlloc(numCols+1, PS_TYPE_F64); 657 950 658 951 // Fill vectors with results and status of results -
trunk/psLib/src/types/psLookupTable.h
r3264 r3271 7 7 * @author Ross Harman, MHPCC 8 8 * 9 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-02-17 19:26:24$9 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-02-17 21:54:09 $ 11 11 * 12 * Copyright 2004 -2005Maui High Performance Computing Center, University of Hawaii12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 13 13 */ 14 14
Note:
See TracChangeset
for help on using the changeset viewer.
