/** @file  psLookupTable.c
*
*  @brief This file defines the structure and functions for table lookups.
*
*  @ingroup types
*
*  @author Ross Harman, MHPCC
*
*  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
*  @date $Date: 2008-04-17 23:43:03 $
*
*  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 "psAbort.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 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;
}

#define PARSE_VALUE_INT_CASE(TYPE, FUNC) \
    case PS_TYPE_##TYPE: { \
        char *end = NULL; \
        ps##TYPE value = FUNC(strValue, &end, 0); \
        if (*end != '\0' && !isspace(*end)) { \
            *status = PS_PARSE_ERROR_VALUE; \
        } \
        vec->data.TYPE[index] = value; \
        return; \
    }

#define PARSE_VALUE_FLOAT_CASE(TYPE, FUNC) \
    case PS_TYPE_##TYPE: { \
        char *end = NULL; \
        ps##TYPE value = FUNC(strValue, &end); \
        if (*end != '\0' && !isspace(*end)) { \
            *status = PS_PARSE_ERROR_VALUE; \
        } \
        vec->data.TYPE[index] = value; \
        return; \
    }


/* 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)
{
    switch (vec->type.type) {
        PARSE_VALUE_INT_CASE(S32, strtol);
        PARSE_VALUE_INT_CASE(S64, strtoll);
        PARSE_VALUE_FLOAT_CASE(F32, strtof);
        PARSE_VALUE_FLOAT_CASE(F64, strtod);
      default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for vector: %x.\n", vec->type.type);
        *status = PS_PARSE_ERROR_TYPE;
    }
    return;
}

static void lookupTableFree(psLookupTable* table)
{
    psAssert(table, "impossible");
    psFree(table->values);
    psFree(table->filename);
    psFree(table->format);
    psFree(table->index);
}

/*****************************************************************************/
/* 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)
{
    PS_ASSERT_STRING_NON_EMPTY(filename,NULL);
    PS_ASSERT_STRING_NON_EMPTY(format,NULL);

    psLookupTable *outTable = 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;
}

psArray *psVectorsReadFromFile(const char *filename,
                               const char *format)
{
    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
    PS_ASSERT_STRING_NON_EMPTY(format, NULL);

    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;

    // 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;
}

bool p_psLookupTableImport(const char *file,
                           unsigned int lineno,
                           const char *func,
                           psLookupTable *table,
                           psArray *vectors,
                           long indexCol)
{
    PS_ASSERT_LOOKUPTABLE_NON_NULL(table, false);
    PS_ASSERT_ARRAY_NON_NULL(vectors, false);
    PS_ASSERT_INT_NONNEGATIVE(indexCol, false);
    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 false;
    }

    psVector *indexVector = vectors->data[indexCol]; // The vector that provides the index
    long numRows = indexVector->n;      // Number of rows in table
    long numCols = vectors->n;          // Number of columns in table

#define CHECK_SORT_CASE(TYPE) \
    case PS_TYPE_##TYPE: \
        for (long i = 1; i < numRows; i++) { \
             if (indexVector->data.TYPE[i] < indexVector->data.TYPE[i - 1]) { \
                 sortNeeded = true; \
                 break; \
             } \
         } \
         break;

    // Check if array is sorted on the index column
    bool sortNeeded  = false;           // Do we need to sort?
    switch (indexVector->type.type) {
        CHECK_SORT_CASE(U8);
        CHECK_SORT_CASE(U16);
        CHECK_SORT_CASE(U32);
        CHECK_SORT_CASE(U64);
        CHECK_SORT_CASE(S8);
        CHECK_SORT_CASE(S16);
        CHECK_SORT_CASE(S32);
        CHECK_SORT_CASE(S64);
        CHECK_SORT_CASE(F32);
        CHECK_SORT_CASE(F64);
      default:
        psAbort("Should never get here: bad type.");
    }

    if (sortNeeded) {
        psVector *sortedIndex = psVectorSortIndex(NULL, indexVector); // Indices for sorting
        psArray *newValues = psArrayAlloc(numCols);
        for (long i = 0; i < numCols; i++) {
            psVector *oldVector = vectors->data[i]; // Vector of interest
            psVector *newVector = psVectorAlloc(numRows, oldVector->type.type);
            newValues->data[i] = newVector;

#define REARRANGE_CASE(TYPE) \
    case PS_TYPE_##TYPE: \
        for (long j = 0; j < numRows; j++) { \
            newVector->data.TYPE[j] = oldVector->data.TYPE[sortedIndex->data.S32[j]]; \
        } \
        break;

            switch (oldVector->type.type) {
                REARRANGE_CASE(U8);
                REARRANGE_CASE(U16);
                REARRANGE_CASE(U32);
                REARRANGE_CASE(U64);
                REARRANGE_CASE(S8);
                REARRANGE_CASE(S16);
                REARRANGE_CASE(S32);
                REARRANGE_CASE(S64);
                REARRANGE_CASE(F32);
                REARRANGE_CASE(F64);
              default:
                psAbort("Should never get here: bad type.");
            }
        }

        psFree(sortedIndex);
        table->values = newValues;
    } else {
        table->values = psMemIncrRefCounter(vectors);
    }

    table->index = psMemIncrRefCounter(table->values->data[indexCol]);
    table->indexCol = indexCol;

#define SET_VALID_CASE(TYPE) \
    case PS_TYPE_##TYPE: \
        *(double*)&table->validFrom = table->index->data.TYPE[0]; \
        *(double*)&table->validTo   = table->index->data.TYPE[numRows - 1]; \
        break;

    switch (table->index->type.type) {
        SET_VALID_CASE(U8);
        SET_VALID_CASE(U16);
        SET_VALID_CASE(U32);
        SET_VALID_CASE(U64);
        SET_VALID_CASE(S8);
        SET_VALID_CASE(S16);
        SET_VALID_CASE(S32);
        SET_VALID_CASE(S64);
        SET_VALID_CASE(F32);
        SET_VALID_CASE(F64);
      default:
        psAbort("Should never get here: bad type.");
    }

    return true;
}

long psLookupTableRead(psLookupTable* table)
{
    PS_ASSERT_LOOKUPTABLE_NON_NULL(table, 0);

    // Read vectors from file specified in table
    psArray *vectors = psVectorsReadFromFile(table->filename, table->format); // Vectors in file
    if (!vectors) {
        psError(PS_ERR_UNKNOWN, false, "Unable to read vectors from %s", table->filename);
        return 0;
    }

    // Import vector into table
    if (!psLookupTableImport(table, vectors, table->indexCol)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to import vectors from %s into lookup table.",
                table->filename);
        psFree(vectors);
        return 0;
    }
    psFree(vectors);                // Drop reference

    return table->index->n;
}

#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)
{
    PS_ASSERT_LOOKUPTABLE_NON_NULL(table, NAN);
    // Ensuring information exists
    PS_ASSERT_VECTOR_NON_NULL(table->index, NAN);
    PS_ASSERT_ARRAY_NON_NULL(table->values, NAN);
    PS_ASSERT_INT_WITHIN_RANGE(column, 0, (int)table->values->n - 1, NAN);
    PS_ASSERT_VECTOR_NON_NULL((psVector*)table->values->data[column], NAN);

    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;

    indexVec = table->index;
    values = table->values;
    numRows = table->index->n;
    numCols = table->values->n;
    valuesVec = (psVector*)values->data[column];

    // 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)
{
    PS_ASSERT_LOOKUPTABLE_NON_NULL(table, NULL);
    // Ensuring information exists
    PS_ASSERT_VECTOR_NON_NULL(table->index, NULL);
    PS_ASSERT_ARRAY_NON_NULL(table->values, NULL);
    long numCols = table->values->n;
    if (numCols == 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No columns in lookup table from which to interpolate.");
        return NULL;
    }

    psVector *outVector = psVectorAlloc(numCols, PS_TYPE_F64);

    // Fill vectors with results and status of results
    // XXX: This algorithm must be changed to something more efficient!
    for (long 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;
}

