/** @file  psLookupTable.c
*
*  @brief This file defines the structure and functions for table lookups.
*
*  @ingroup types
*
*  @author Ross Harman, MHPCC
*
*  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
*  @date $Date: 2007-02-08 21:57:02 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii
*/

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <stdio.h>
#include <string.h>
#include <ctype.h>
#undef __STRICT_ANSI__
#include <stdlib.h>
#define __STRICT_ANSI__
#include <limits.h>
#include <math.h>
#include <inttypes.h>

#include "psMemory.h"
#include "psString.h"
#include "psError.h"
#include "psLookupTable.h"

#include "psAssert.h"

/******************************************************************************/
/*  DEFINE STATEMENTS                                                         */
/******************************************************************************/

/** Maximum size of a string */
#define MAX_STRING_LENGTH   256
#define ARRAY_STRIDE        16
#define INITIAL_NUM         10          // Initial number of elements
/******************************************************************************/
/*  TYPE DEFINITIONS                                                          */
/******************************************************************************/

// None

/*****************************************************************************/
/*  GLOBAL VARIABLES                                                         */
/*****************************************************************************/

// None

/*****************************************************************************/
/*  FILE STATIC VARIABLES                                                    */
/*****************************************************************************/

// None

/*****************************************************************************/
/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
/*****************************************************************************/

static bool ignoreLine(char *inString);
static char *cleanString(char *inString, int sLen);
static char* getToken(char **inString, char *delimiter, psParseErrorType *status);
static psS32 parseS32(char *inString, psParseErrorType *status);
static psS64 parseS64(char *inString, psParseErrorType *status);
static psF32 parseF32(char *inString, psParseErrorType *status);
static psF64 parseF64(char *inString, psParseErrorType *status);
static void parseValue(psVector *vec, psU64 index, char* strValue, psParseErrorType *status);
static void lookupTableFree(psLookupTable* table);

/* Determines if a line is blank (whitespace only) or a commentline. It returns true if so.
   The input string must be null terminated. */
static bool ignoreLine(char *inString)
{
    while(*inString!='\0' && *inString!='#') {
        if(!isspace(*inString)) {
            return false;
        }
        inString++;
    }
    return true;
}


/* Removes leading and trailing whitespace and # characters from a string. The
   cleaned string is a new null terminated copy of the original input string. */
static char *cleanString(char *inString,
                         int sLen)
{
    char *ptrB = NULL;
    char *ptrE = NULL;
    char *cleaned = NULL;

    ptrB = inString;

    // Skip over leading # or whitespace
    //    while (isspace(*ptrB) || *ptrB=='#') {
    //        ptrB++;
    //    }

    // Skip over trailing whitespace, null terminators, and # characters
    ptrE = inString + sLen;
    while(isspace(*ptrE) || *ptrE=='\0' || *ptrE=='#') {
        ptrE--;
    }

    // Length, sLen, does not include '\0'
    sLen = ptrE - ptrB + 1;

    // Adds '\0' to end of string and +1 to sLen
    cleaned = psStringNCopy(ptrB, sLen);

    return cleaned;
}


/* Returns cleaned token based on delimiter, but not including delimiter. Also changes
   the pointer location the beginning of the string. Tokens are newly allocated null
    terminated strings. */
static char* getToken(char **inString,
                      char *delimiter,
                      psParseErrorType *status)
{
    char *cleanToken = NULL;
    int sLen = 0;

    // Skip over leading whitespace
    while(isspace(**inString)) {
        (*inString)++;
    }

    // Length of token, not including delimiter
    sLen = strcspn(*inString, delimiter);
    if(sLen) {

        // Create new, cleaned, and null terminated token
        cleanToken = cleanString(*inString, sLen);

        // Move to end of token
        (*inString) += sLen;
    }
    /*    else if(**inString!='\0' && sLen==0) {
            *status = PS_PARSE_ERROR_GENERAL;
        }
    */
    return cleanToken;
}

/* Returns single parsed value as a psS32. The input string must be cleaned and
   null terminated. */
static psS32 parseS32(char *inString,
                      psParseErrorType *status)
{
    char *end = NULL;
    psS32 value = 0.0;


    value = (psS32)strtol(inString, &end, 0);
    if(*end != '\0' && !isspace(*end)) {
        *status = PS_PARSE_ERROR_VALUE;
    }
    /*    else if(inString==end) {
            *status = PS_PARSE_ERROR_VALUE;
        }
    */
    return value;
}

/* Returns single parsed value as a psS64. The input string must be cleaned and
   null terminated. */
static psS64 parseS64(char *inString,
                      psParseErrorType *status)
{
    char *end = NULL;
    psS64 value = 0.0;

    value = (psS64)strtoll(inString, &end, 0);
    if(*end != '\0' && !isspace(*end)) {
        *status = PS_PARSE_ERROR_VALUE;
    }
    /*    else if(inString==end) {
            *status = PS_PARSE_ERROR_VALUE;
        }
    */
    return value;
}

/* Returns single parsed value as a psF32. The input string must be cleaned and
   null terminated. */
static psF32 parseF32(char *inString,
                      psParseErrorType *status)
{
    char *end = NULL;
    psF32 value = 0.0;

    value = (psF32)strtof(inString, &end);
    if(*end != '\0' && !isspace(*end)) {
        *status = PS_PARSE_ERROR_VALUE;
    }
    /*    else if(inString==end) {
            *status = PS_PARSE_ERROR_VALUE;
        }
    */
    return value;
}

/* Returns single parsed value as a psF64. The input string must be cleaned and
   null terminated. */
static psF64 parseF64(char *inString,
                      psParseErrorType *status)
{
    char *end = NULL;
    psF64 value = 0.0;

    value = (psF64)strtod(inString, &end);
    if(*end != '\0' && !isspace(*end)) {
        *status = PS_PARSE_ERROR_VALUE;
    }
    /*    else if(inString==end) {
            *status = PS_PARSE_ERROR_VALUE;
        }
    */
    return value;
}

/* Returns single parsed value as a double precision number. The input string must be
   cleaned and null terminated. */
static void parseValue(psVector *vec,
                       psU64 index,
                       char* strValue,
                       psParseErrorType *status)
{
    psElemType type;


    /*    if(vec == NULL) {
            *status = 1;
            return;
        }
    */
    type = vec->type.type;

    switch(type) {
    case PS_TYPE_S32:
        vec->data.S32[index] = parseS32(strValue, status);
        break;
    case PS_TYPE_S64:
        vec->data.S64[index] = parseS64(strValue, status);
        break;
    case PS_TYPE_F32:
        vec->data.F32[index] = parseF32(strValue, status);
        break;
    case PS_TYPE_F64:
    default:
        vec->data.F64[index] = parseF64(strValue, status);
        break;
        //    default:
        //        *status = PS_PARSE_ERROR_TYPE;
    }

    return;
}

static void lookupTableFree(psLookupTable* table)
{
    psFree(table->values);
    psFree((void *)table->filename);
    psFree((void *)table->format);
}

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
/*****************************************************************************/

bool psMemCheckLookupTable(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)lookupTableFree );
}


psLookupTable* p_psLookupTableAlloc(const char *file,
                                    unsigned int lineno,
                                    const char *func,
                                    const char *filename,
                                    const char *format,
                                    long indexCol)
{
    // Can't read table if you don't know its name
    PS_ASSERT_PTR_NON_NULL(filename,NULL);

    // Can't read table if you don't know its format
    PS_ASSERT_PTR_NON_NULL(format,NULL);

    psLookupTable *outTable = NULL;

    // Allocate lookup table
    outTable = (psLookupTable*)p_psAlloc(file, lineno, func, sizeof(psLookupTable));

    // Set deallocator
    psMemSetDeallocator(outTable, (psFreeFunc)lookupTableFree);

    // Allocate and set file name and format strings
    outTable->filename = psStringCopy(filename);
    outTable->format = psStringCopy(format);

    // Valid ranges. Automatically set by table read if both zero.
    *(double *)&outTable->validFrom = 0;
    *(double *)&outTable->validTo = 0;
    outTable->indexCol = indexCol;

    // Vector of independent index values. Filled by table read.
    outTable->index = NULL;

    // Array of dependent table values corresponding to index values. Filled by table read.
    outTable->values = NULL;

    return outTable;
}

#define UPDATE_VALID_TO_FROM(TABLE)                                             \
switch (TABLE->index->type.type) {                                              \
case PS_TYPE_U8:                                                                \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.U8[0];                         \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.U8[TABLE->index->n-1];         \
    break;                                                                      \
case PS_TYPE_S8:                                                                \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.S8[0];                         \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.S8[TABLE->index->n-1];         \
    break;                                                                      \
case PS_TYPE_U16:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.U16[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.U16[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_S16:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.S16[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.S16[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_U32:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.U32[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.U32[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_S32:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.S32[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.S32[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_U64:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.U64[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.U64[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_S64:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.S64[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.S64[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_F32:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.F32[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.F32[TABLE->index->n-1];        \
    break;                                                                      \
case PS_TYPE_F64:                                                               \
    *(double *)&TABLE->validFrom = (psF64)TABLE->index->data.F64[0];                        \
    *(double *)&TABLE->validTo   = (psF64)TABLE->index->data.F64[TABLE->index->n-1];        \
    break;                                                                      \
default:                                                                        \
    *(double *)&TABLE->validFrom = (psF64)0;                                                \
    *(double *)&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];   \
    if(INDEX_OUT == 0) {                                                                   \
        validFrom = ((psVector*)VEC_OUT)->data.S32[INDEX_OUT];                             \
    }                                                                                      \
    if(INDEX_OUT == ((psVector*)VEC_OUT)->n-1 ) {                                          \
        validTo = ((psVector*)VEC_OUT)->data.S32[INDEX_OUT];                               \
    }                                                                                      \
    break;                                                                                 \
case PS_TYPE_S64:                                                                          \
    ((psVector*)VEC_OUT)->data.S64[INDEX_OUT] = ((psVector*)VEC_IN)->data.S64[INDEX_IN];   \
    if(INDEX_OUT == 0) {                                                                   \
        validFrom = ((psVector*)VEC_OUT)->data.S64[INDEX_OUT];                             \
    }                                                                                      \
    if(INDEX_OUT == ((psVector*)VEC_OUT)->n-1 ) {                                          \
        validTo = ((psVector*)VEC_OUT)->data.S64[INDEX_OUT];                               \
    }                                                                                      \
    break;                                                                                 \
case PS_TYPE_F32:                                                                          \
    ((psVector*)VEC_OUT)->data.F32[INDEX_OUT] = ((psVector*)VEC_IN)->data.F32[INDEX_IN];   \
    if(INDEX_OUT == 0) {                                                                   \
        validFrom = ((psVector*)VEC_OUT)->data.F32[INDEX_OUT];                             \
    }                                                                                      \
    if(INDEX_OUT == ((psVector*)VEC_OUT)->n-1 ) {                                          \
        validTo = ((psVector*)VEC_OUT)->data.F32[INDEX_OUT];                               \
    }                                                                                      \
    break;                                                                                 \
case PS_TYPE_F64:                                                                          \
    ((psVector*)VEC_OUT)->data.F64[INDEX_OUT] = ((psVector*)VEC_IN)->data.F64[INDEX_IN];   \
    if(INDEX_OUT == 0) {                                                                   \
        validFrom = ((psVector*)VEC_OUT)->data.F64[INDEX_OUT];                             \
    }                                                                                      \
    if(INDEX_OUT == ((psVector*)VEC_OUT)->n-1 ) {                                          \
        validTo = ((psVector*)VEC_OUT)->data.F64[INDEX_OUT];                               \
    }                                                                                      \
    break;                                                                                 \
default:                                                                                   \
    break;                                                                                 \
}

psArray *psVectorsReadFromFile(const char *filename,
                               const char *format)
{
    psArray*          outputArray = NULL;
    psVector*         colVector   = NULL;
    char*             strValue    = NULL;
    char*             strNum      = NULL;
    char*             line        = NULL;
    char*             linePtr     = NULL;
    int               numCols     = 0;
    int               numRows     = 0;
    FILE*             fp          = NULL;
    const char*       tempFormat  = NULL;
    psParseErrorType  parseStatus = PS_PARSE_SUCCESS;

    // Check for nul file name
    PS_ASSERT_PTR_NON_NULL(filename,NULL);

    // Parse format string for valid string
    PS_ASSERT_PTR_NON_NULL(format,NULL);

    // Create temp pointer which can then be used several times
    tempFormat = format;

    // Create output array and set array elements to zero
    outputArray = psArrayAllocEmpty(INITIAL_NUM);

    // Parse the format string to determine how many vectors
    // and whether the format string is valid
    while((strValue=getToken((char**)&tempFormat," \t",&parseStatus))) {

        // Check for %d format sub string
        if(strcmp(strValue,"\%d") == 0 ) {
            numCols++;
            colVector = psVectorAlloc(1,PS_TYPE_S32);
            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
            psFree(colVector);
        } else if (strcmp(strValue,"\%ld") == 0) {
            numCols++;
            colVector = psVectorAlloc(1,PS_TYPE_S64);
            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
            psFree(colVector);
        } else if (strcmp(strValue,"\%f") == 0) {
            numCols++;
            colVector = psVectorAlloc(1,PS_TYPE_F32);
            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
            psFree(colVector);
        } else if (strcmp(strValue,"\%lf") == 0) {
            numCols++;
            colVector = psVectorAlloc(1,PS_TYPE_F64);
            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
            psFree(colVector);
        } else if (strstr(strValue,"\%*") != 0) {
            // Don't increase number of columns
        } else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    "Invalid format specifier");
            psFree(strValue);
            numCols = 0;
            break;
        }
        psFree(strValue);
    }

    // If the format string was parsed successfully and return numCols the
    // prepare to open file and read values
    if(numCols > 0) {

        // Open specified file
        if((fp=fopen(filename, "r")) == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Failed to open file %s."),
                    filename);
            psFree(outputArray);
            return NULL;
        } else {
            // Initialize array index
            int arrayIndex = 0;

            // Create reusable line for continuous read
            line = (char*)psAlloc(MAX_STRING_LENGTH*sizeof(char));

            // Loop through file to get numRows, numCols, and column data types
            while((fgets(line, MAX_STRING_LENGTH, fp) != NULL) &&
                    (parseStatus == PS_PARSE_SUCCESS))   {

                // Copy pointer to line for parsing
                linePtr = line;

                // If line is not a comment or blank, then extract data
                if(!ignoreLine(linePtr)) {
                    numRows++;

                    // Copy format pointer for parsing
                    tempFormat = format;
                    arrayIndex = 0;
                    parseStatus = PS_PARSE_SUCCESS;

                    // Loop through format and line strings to get values in text table file
                    while((strValue=getToken((char**)&tempFormat," \t",&parseStatus))
                            && (strNum=getToken((char**)&linePtr," \t",&parseStatus)) ) {
                        // Set column vector
                        colVector = outputArray->data[arrayIndex];

                        // Set column entries based on format string defining the type
                        if(strcmp(strValue,"\%d") == 0 ) {
                            colVector = psVectorRecycle(colVector, numRows,
                                                        colVector->type.type);
                            parseValue(colVector,numRows-1,strNum,&parseStatus);
                            arrayIndex++;
                        } else if (strcmp(strValue,"\%ld") == 0) {
                            colVector = psVectorRecycle(colVector, numRows,
                                                        colVector->type.type);
                            parseValue(colVector,numRows-1,strNum,&parseStatus);
                            arrayIndex++;
                        } else if (strcmp(strValue,"\%f") == 0) {
                            colVector = psVectorRecycle(colVector, numRows,
                                                        colVector->type.type);
                            parseValue(colVector,numRows-1,strNum,&parseStatus);
                            arrayIndex++;
                        } else if (strcmp(strValue,"\%lf") == 0) {
                            colVector = psVectorRecycle(colVector, numRows,
                                                        colVector->type.type);
                            parseValue(colVector,numRows-1,strNum,&parseStatus);
                            arrayIndex++;
                        } else if (strstr(strValue,"\%*") != 0) {
                            // Don't increase number of columns
                        }
                        psFree(strValue);
                        psFree(strNum);

                        // If the file line was not parsed successful report
                        // error and return NULL
                        if(parseStatus != PS_PARSE_SUCCESS) {
                            psError(PS_ERR_UNKNOWN, true,
                                    "Parsing text file failed.");
                            fclose(fp);
                            psFree(outputArray);
                            psFree(line);
                            return NULL;
                        }
                    }
                    if (strValue != NULL && strNum == NULL) {
                        psError(PS_ERR_UNKNOWN, true,
                                "Parsing text file failed - missing table value(s).");
                        fclose(fp);
                        psFree(outputArray);
                        psFree(line);
                        psFree(strValue);
                        return NULL;
                    }
                }  // ignore line
            }

            //Return NULL for an empty table
            if (numRows == 0) {
                psError(PS_ERR_UNKNOWN, true,
                        "Parsing text file failed - input table is empty.");
                fclose(fp);
                psFree(outputArray);
                psFree(line);
                return NULL;
            }

            // Read on the lines in the file - close file pointer
            fclose(fp);
            psFree(line);
        }
    } else {
        // Format string parse error detected
        psError(PS_ERR_UNKNOWN, true,
                "Format string was not parsed sucessfully");
        psFree(outputArray);
        return NULL;
    }

    // Return populated array
    return outputArray;
}

psLookupTable *p_psLookupTableImport(const char *file,
                                     unsigned int lineno,
                                     const char *func,
                                     psLookupTable *table,
                                     const psArray *vectors,
                                     long indexCol)
{
    psLookupTable* outputTable = NULL;
    bool         sortNeeded  = false;
    psF64          validTo     = 0;
    psF64          validFrom   = 0;

    // Check for NULL input table
    // XXX table/outputtable should be allocated if NULL
    PS_ASSERT_PTR_NON_NULL(table,NULL);

    // Check for NULL vectors
    PS_ASSERT_PTR_NON_NULL(vectors,NULL);

    // Check for invalid index column
    if(indexCol < 0 ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "Index column cannot be less than zero");
        return NULL;
    }

    if (indexCol >= vectors->n) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Index column, %ld, is larger"
                " than number of columns, %ld.", indexCol, vectors->n);
        return NULL;
    }

    psVector* indexColumnVector = vectors->data[indexCol];
    psS32  numRows = indexColumnVector->n;
    psS32  numCols = vectors->n;

    // Check if array is sorted on the index column
    psVector* sortedIndex = psVectorAlloc(numRows,PS_TYPE_U32);
    sortedIndex = psVectorSortIndex(sortedIndex,indexColumnVector);
    for(psS32 i = 0; i < numRows; i++ ) {
        if(sortedIndex->data.U32[i] != i) {
            sortNeeded = true;
            break;
        }
    }

    // Check if it is necessary to sort value vectors
    if(sortNeeded) {
        // Allocate new array to contain sorted vectors
        psArray* newValueArray = psArrayAlloc(numCols);
        for(psS32 j = 0; j < numCols; j++) {
            psS32 type = ((psVector*)(vectors->data[j]))->type.type;
            newValueArray->data[j] = psVectorAlloc(numRows,type);
        }

        // Populate new array/vectors with sorted values
        for(psS32 i = 0; i < numRows; i++) {
            // Populate new index vector
            psU32 sortIndex = sortedIndex->data.U32[i];
            // For every column populate new value vectors
            for(psS32 j=0; j < numCols; j++) {
                COPY_VECTOR_VALUES(newValueArray->data[j],i,vectors->data[j], sortIndex)
            }
        }
        // Assign new vector value to table
        table->index = newValueArray->data[indexCol];
        // Assign new value vectors to table array
        table->values = newValueArray;
        // Assign indexCol
        table->indexCol = indexCol;
        // Assign validTo and validFrom values
        UPDATE_VALID_TO_FROM(table)

        // Set return value
        outputTable = table;

    } else {
        // Index vector is already sorted
        // Assign array vector specified by indexCol to table index vector
        // and increment memory ref
        table->index = vectors->data[indexCol];
        table->values = (psArray*)vectors;
        psMemIncrRefCounter((psArray*)vectors);

        // Assign indexCol
        table->indexCol = indexCol;
        // Assign validTo and validFrom values
        UPDATE_VALID_TO_FROM(table);

        // Set return value
        outputTable = table;
    }

    // Free sort vector
    psFree(sortedIndex);

    return outputTable;
}

long psLookupTableRead(psLookupTable* table)
{
    long            numRows  = 0;
    psArray*        vectors  = NULL;
    psLookupTable*  outTable = NULL;

    // Check if the input table is NULL then return 0
    PS_ASSERT_PTR_NON_NULL(table,0);

    // Read vectors from file specified in table
    vectors = psVectorsReadFromFile(table->filename, table->format);

    // Check for valid array after reading from file
    if(vectors != NULL) {

        // Import vector into table
        outTable = psLookupTableImport(table, vectors, table->indexCol);

        psFree(vectors);

        // Update the number of rows read if outTable is not NULL
        if(outTable != NULL) {
            numRows = table->index->n;
        }
    }

    // Return the number of lines read
    return numRows;
}

#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_S32:                                                       \
    if( (psS32)index < TABLE->index->data.S32[0] ) {                    \
        return NAN;                                                     \
    }                                                                   \
    if( (psS32)index > TABLE->index->data.S32[numRows-1] ) {            \
        return NAN;                                                     \
    }                                                                   \
    break;                                                              \
case PS_TYPE_S64:                                                       \
    if( (psS64)index < TABLE->index->data.S64[0] ) {                    \
        return NAN;                                                     \
    }                                                                   \
    if( (psS64)index > TABLE->index->data.S64[numRows-1] ) {            \
        return NAN;                                                     \
    }                                                                   \
    break;                                                              \
case PS_TYPE_F32:                                                       \
    if( (psF32)index < TABLE->index->data.F32[0] ) {                    \
        return NAN;                                                     \
    }                                                                   \
    if( (psF32)index > TABLE->index->data.F32[numRows-1] ) {            \
        return NAN;                                                     \
    }                                                                   \
    break;                                                              \
case PS_TYPE_F64:                                                       \
    if( index < TABLE->index->data.F64[0] ) {                           \
        return NAN;                                                     \
    }                                                                   \
    if( index > TABLE->index->data.F64[numRows-1] ) {                   \
        return NAN;                                                     \
    }                                                                   \
    break;                                                              \
default:                                                                \
    return NAN;                                                         \
    break;                                                              \
}

double psLookupTableInterpolate(const psLookupTable *table,
                                double index,
                                long column)
{
    psU64 hiIdx = 0;
    psU64 loIdx = 0;
    long  numRows = 0;
    long numCols = 0;
    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;

    // Check for NULL table
    PS_ASSERT_PTR_NON_NULL(table,NAN);

    indexVec = table->index;
    values = table->values;
    numRows = table->index->n;
    numCols = table->values->n;
    PS_ASSERT_PTR_NON_NULL(indexVec,NAN);
    PS_ASSERT_PTR_NON_NULL(values,NAN);
    PS_ASSERT_INT_UNEQUAL(numRows, 0,NAN);
    PS_ASSERT_INT_UNEQUAL(numCols, 0,NAN);
    PS_ASSERT_INT_WITHIN_RANGE(column, 0, (int)(numCols-1), NAN);

    valuesVec = (psVector*)values->data[column];
    PS_ASSERT_PTR_NON_NULL(indexVec,NAN);

    // 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++;
        /*  XXX:  following is unreachable.
                if(hiIdx >= numRows) {
                    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                            _("High index too big, %" PRIu64 "."), hiIdx);
                    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,
                _("Low index too small, %" PRIu64 "."), loIdx);
        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,
                _("Divide by zero error during interpolation."));
        return NAN;
    } else {
        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;
    }

    return out;
}

psVector* psLookupTableInterpolateAll(const psLookupTable *table,
                                      double index)
{
    long numCols = 0;
    psVector *outVector = NULL;

    // Error checks
    PS_ASSERT_PTR_NON_NULL(table,NULL);
    PS_ASSERT_PTR_NON_NULL(table->values,NULL);
    numCols = table->values->n;
    PS_ASSERT_INT_UNEQUAL(numCols, 0,NULL);

    // Create output vector
    outVector = psVectorAlloc(numCols, PS_TYPE_F64);

    // Fill vectors with results and status of results
    for(psU64 i=0; i<numCols; i++) {
        outVector->data.F64[i] = psLookupTableInterpolate(table, index, i);
        if(isnan(outVector->data.F64[i])) {
            // Free allocated vector
            psFree(outVector);
            outVector = NULL;
            // Break out of loop since error detected
            break;
        }
    }

    return outVector;
}

