/** @file  psVector.c
*
*  @brief Contains support for basic vector types
*
*  This file defines the basic type for a vector struct and functions useful
*  in manupulating vectors.
*
*  @author Ross Harman, MHPCC
*  @author Robert DeSonia, MHPCC
*
*  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-10-04 01:31:40 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#include <string.h>                        // for memcpy
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <unistd.h>

#include "psMemory.h"
#include "psError.h"
#include "psVector.h"
#include "psLogMsg.h"
#include "psCompare.h"

#include "psErrorText.h"

typedef struct
{
    union {
        psU8* U8;                      ///< Unsigned 8-bit integer data.
        psU16* U16;                    ///< Unsigned 16-bit integer data.
        psU32* U32;                    ///< Unsigned 32-bit integer data.
        psU64* U64;                    ///< Unsigned 64-bit integer data.
        psS8* S8;                      ///< Signed 8-bit integer data.
        psS16* S16;                    ///< Signed 16-bit integer data.
        psS32* S32;                    ///< Signed 32-bit integer data.
        psS64* S64;                    ///< Signed 64-bit integer data.
        psF32* F32;                    ///< Single-precision float data.
        psF64* F64;                    ///< Double-precision float data.
        psC32* C32;                    ///< Single-precision complex data.
        psC64* C64;                    ///< Double-precision complex data.
    } data;
    psU32 index;
}
indexedVector;

static void vectorFree(psVector* psVec);

static void vectorFree(psVector* psVec)
{
    if (psVec == NULL) {
        return;
    }

    psFree(psVec->data.U8);
}

// FUNCTION IMPLEMENTATION - PUBLIC
bool psMemCheckVector(psPtr ptr)
{
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)vectorFree );
}


psVector* psVectorAlloc(long nalloc,
                        psElemType type)
{
    psVector* psVec = NULL;
    psS32 elementSize = 0;

    elementSize = PSELEMTYPE_SIZEOF(type);
    if (elementSize < 1) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psVector_UNSUPPORTED_TYPE, type);
        return NULL;
    }


    // Create vector struct
    psVec = (psVector* ) psAlloc(sizeof(psVector));
    psMemSetDeallocator(psVec, (psFreeFunc) vectorFree);

    psVec->type.dimen = PS_DIMEN_VECTOR;
    psVec->type.type = type;
    *(long*)&psVec->nalloc = nalloc;
    psVec->n = nalloc;

    // Create vector data array
    psVec->data.U8 = psAlloc(nalloc * elementSize);

    return psVec;
}

psVector* psVectorRealloc(psVector* vector,
                          long nalloc)
{
    psS32 elementSize = 0;
    psElemType elemType;

    if (vector == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psVector_REALLOC_NULL);
        return NULL;
    } else if (vector->nalloc != nalloc) {     // No need to realloc to same size
        elemType = vector->type.type;
        elementSize = PSELEMTYPE_SIZEOF(elemType);
        if (nalloc < vector->n) {
            vector->n = nalloc;
        }
        // Realloc after decrementation to avoid accessing freed array elements
        vector->data.U8 = psRealloc(vector->data.U8, nalloc * elementSize);
        *(long*)&vector->nalloc = nalloc;
    }

    return vector;
}

psVector* psVectorRecycle(psVector* vector,
                          long nalloc,
                          psElemType type)
{
    psS32 byteSize;

    if (vector == NULL) {
        return psVectorAlloc(nalloc, type);
    }

    if (vector->type.dimen !=  PS_DIMEN_VECTOR &&
            vector->type.dimen !=  PS_DIMEN_TRANSV) {
        psFree(vector);
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psVector_NOT_A_VECTOR);
        return NULL;
    }

    byteSize = nalloc * PSELEMTYPE_SIZEOF(type);

    // need to increase data buffer?
    if (byteSize > vector->nalloc*PSELEMTYPE_SIZEOF(vector->type.type)) {
        vector->data.U8 = psRealloc(vector->data.U8, byteSize);
        *(long*)&vector->nalloc = nalloc;
    }

    vector->type.dimen = PS_DIMEN_VECTOR;
    vector->type.type = type;
    vector->n = nalloc;
    return vector;
}

psVector *psVectorExtend(psVector *vector,
                         long delta,
                         long nExtend)
{
    // can't handle a NULL vector (don't know the data type to allocate)
    if (vector == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psVector_NULL);
        return NULL;
    }

    // requirement on delta; if < 1, set to 10.
    if (delta < 1) {
        delta = 10;
    }

    // adjust the allocated size, if needed. (if nExtend < 1, this is will never happen)
    if (nExtend > 0) {
        unsigned int minAlloc = vector->n + nExtend + nExtend;
        if (vector->nalloc < minAlloc) {
            unsigned int nAlloc = delta + vector->nalloc;
            // make sure the delta is large enough hold twice the extended length.
            if (nAlloc < minAlloc) {
                nAlloc = minAlloc;
            }
            vector = psVectorRealloc(vector, nAlloc);
        }
    } else if (nExtend < -vector->n) {
        // For the case of a negative nExtend, need to check that we are not decreasing
        // vector beyond its own length (i.e., creating a negative length).
        nExtend = -vector->n;
    }

    // increment the length by the value specified
    vector->n += nExtend;

    return vector;
}

psVector* psVectorCopy(psVector* output,
                       const psVector* input,
                       psElemType type)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psVector_SORT_NULL);
        psFree(output);
        return NULL;
    }

    psS32 nElements = input->n;

    output = psVectorRecycle(output, nElements, type);

    #define PSVECTOR_COPY(INTYPE,OUTTYPE) { \
        ps##INTYPE *inVec = input->data.INTYPE; \
        ps##OUTTYPE *outVec = output->data.OUTTYPE; \
        for (psS32 col=0;col<nElements;col++) { \
            *(outVec++) = *(inVec++); \
        } \
    }

    #define PSVECTOR_COPY_CASE(OUTTYPE) \
case PS_TYPE_##OUTTYPE: { \
        switch (input->type.type) { \
        case PS_TYPE_S8: \
            PSVECTOR_COPY(S8,OUTTYPE); \
            break; \
        case PS_TYPE_S16: \
            PSVECTOR_COPY(S16,OUTTYPE); \
            break; \
        case PS_TYPE_S32: \
            PSVECTOR_COPY(S32,OUTTYPE); \
            break; \
        case PS_TYPE_S64: \
            PSVECTOR_COPY(S64,OUTTYPE); \
            break; \
        case PS_TYPE_U8: \
            PSVECTOR_COPY(U8,OUTTYPE); \
            break; \
        case PS_TYPE_U16: \
            PSVECTOR_COPY(U16,OUTTYPE); \
            break; \
        case PS_TYPE_U32: \
            PSVECTOR_COPY(U32,OUTTYPE); \
            break; \
        case PS_TYPE_U64: \
            PSVECTOR_COPY(U64,OUTTYPE); \
            break; \
        case PS_TYPE_F32: \
            PSVECTOR_COPY(F32,OUTTYPE); \
            break; \
        case PS_TYPE_F64: \
            PSVECTOR_COPY(F64,OUTTYPE); \
            break; \
        case PS_TYPE_C32: \
            PSVECTOR_COPY(C32,OUTTYPE); \
            break; \
        case PS_TYPE_C64: \
            PSVECTOR_COPY(C64,OUTTYPE); \
            break; \
        default: { \
                char* typeStr; \
                PS_TYPE_NAME(typeStr,type); \
                psError(PS_ERR_BAD_PARAMETER_TYPE, true, \
                        PS_ERRORTEXT_psVector_UNSUPPORTED_TYPE, \
                        typeStr); \
                psFree(output); \
            } \
        } \
        break; \
    }

    switch (type) {
        PSVECTOR_COPY_CASE(S8);
        PSVECTOR_COPY_CASE(S16);
        PSVECTOR_COPY_CASE(S32);
        PSVECTOR_COPY_CASE(S64);
        PSVECTOR_COPY_CASE(U8);
        PSVECTOR_COPY_CASE(U16);
        PSVECTOR_COPY_CASE(U32);
        PSVECTOR_COPY_CASE(U64);
        PSVECTOR_COPY_CASE(F32);
        PSVECTOR_COPY_CASE(F64);
        PSVECTOR_COPY_CASE(C32);
        PSVECTOR_COPY_CASE(C64);
    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psVector_UNSUPPORTED_TYPE,
                    typeStr);
            psFree(output);

            break;
        }
    }
    return output;


}

psVector* psVectorSort(psVector* outVector,
                       const psVector* inVector)
{
    psS32 N = 0;
    psS32 elSize = 0;
    psPtr inVec = NULL;
    psPtr outVec = NULL;
    psElemType inType = 0;

    if (inVector == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psVector_SORT_NULL);
        psFree(outVector);
        return NULL;
    }

    inType = inVector->type.type;
    N = inVector->n;
    inVec = (psPtr)inVector->data.U8;
    elSize = PSELEMTYPE_SIZEOF(inType);

    if (outVector == NULL) {
        outVector = psVectorAlloc(N, inType);
    }

    // check to see if output vector needs to be resized/retyped
    if ( (N > outVector->nalloc) ||
            (inType != outVector->type.type) ) {
        // reshape the output vector to match the input vector's size/type.
        outVector = psVectorRecycle(outVector,N,inType);
    }
    outVector->n = N;
    outVec = outVector->data.U8;

    if (N == 0) {
        // no need to sort anything, as there are no elements in input vector.
        return outVector;
    }

    // Copy input vector values into output vector if not in-place sorting
    if (inVector != outVector) {
        memcpy(outVec, inVec, elSize * N);
    }

    // Sort output vector
    switch (inType) {
    case PS_TYPE_U8:
        qsort(outVec, N, elSize, psCompareU8);
        break;
    case PS_TYPE_U16:
        qsort(outVec, N, elSize, psCompareU16);
        break;
    case PS_TYPE_U32:
        qsort(outVec, N, elSize, psCompareU32);
        break;
    case PS_TYPE_U64:
        qsort(outVec, N, elSize, psCompareU64);
        break;
    case PS_TYPE_S8:
        qsort(outVec, N, elSize, psCompareS8);
        break;
    case PS_TYPE_S16:
        qsort(outVec, N, elSize, psCompareS16);
        break;
    case PS_TYPE_S32:
        qsort(outVec, N, elSize, psCompareS32);
        break;
    case PS_TYPE_S64:
        qsort(outVec, N, elSize, psCompareS64);
        break;
    case PS_TYPE_F32:
        qsort(outVec, N, elSize, psCompareF32);
        break;
    case PS_TYPE_F64:
        qsort(outVec, N, elSize, psCompareF64);
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psVector_UNSUPPORTED_TYPE,
                inType);
        psFree(outVector);
        return NULL;
    }

    return outVector;
}

psVector* psVectorSortIndex(psVector* outVector,
                            const psVector* inVector)
{
    psS32 N = 0;
    psElemType inType = 0;

    if (inVector == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psVector_SORT_NULL);
        psFree(outVector);
        return NULL;
    }

    inType = inVector->type.type;
    N = inVector->n;

    if (N == 0) {
        // no need to sort anything, as there are no elements in input vector.
        return outVector;
    }

    // ok, let's create a temporary indexed vector
    indexedVector* idxVector = psAlloc(sizeof(indexedVector)*N);
    int elSize = PSELEMTYPE_SIZEOF(inType);
    for (int i = 0; i < N; i++) {
        idxVector[i].data.U8 = inVector->data.U8+i*elSize;
        idxVector[i].index = i;
    }

    // Sort indexed vector
    // n.b., since first element in indexedVector is a pointer to the data,
    // we can use the 'Ptr' version of the standard compare functions
    switch (inType) {
    case PS_TYPE_U8:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareU8Ptr);
        break;
    case PS_TYPE_U16:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareU16Ptr);
        break;
    case PS_TYPE_U32:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareU32Ptr);
        break;
    case PS_TYPE_U64:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareU64Ptr);
        break;
    case PS_TYPE_S8:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareS8Ptr);
        break;
    case PS_TYPE_S16:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareS16Ptr);
        break;
    case PS_TYPE_S32:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareS32Ptr);
        break;
    case PS_TYPE_S64:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareS64Ptr);
        break;
    case PS_TYPE_F32:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareF32Ptr);
        break;
    case PS_TYPE_F64:
        qsort(idxVector, N, sizeof(indexedVector), (psCompareFcn)psCompareF64Ptr);
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psVector_UNSUPPORTED_TYPE,
                inType);
        psFree(idxVector);
        psFree(outVector);
        return NULL;
    }

    // extract the indices to the output vector
    outVector = psVectorRecycle(outVector, N, PS_TYPE_U32);
    psU32* outData = outVector->data.U32;
    for (int i = 0; i < N; i++) {
        outData[i] = idxVector[i].index;
    }

    // Free temp memory
    psFree(idxVector);

    return outVector;
}

char* psVectorToString(psVector* vector,
                       int maxLength)
{

    if (maxLength < 5) {
        return NULL;
    }

    char* str = psAlloc(sizeof(char)*maxLength+1);

    if (vector == NULL) {
        snprintf(str,maxLength, "NULL");
        return str;
    }

    int size = vector->n;

    if (size == 0) {
        snprintf(str,maxLength, "[]");
        return str;
    }

    char* tempStr = psAlloc(sizeof(char)*maxLength+1);
    *str = '\0';
    bool full = false;

    #define APPEND_ELEMENTS_CASE(TYPE, NATIVE_TYPE, FORMAT) \
case PS_TYPE_##TYPE: \
    for (lcv=0; lcv < size && ! full; lcv++) { \
        snprintf(tempStr, maxLength, "%s" FORMAT, prefix, (NATIVE_TYPE) (vector->data.TYPE[lcv])); \
        strncat(str,tempStr,maxLength); \
        full = (strlen(str) > maxLength-2); \
        prefix = ","; \
    } \
    break;

    #define APPEND_ELEMENTS_CASE_COMPLEX(TYPE,CREAL,CIMAG) \
case PS_TYPE_##TYPE: \
    for (lcv=0; lcv < size && ! full; lcv++) { \
        snprintf(tempStr, maxLength, "%s%g%+gi", prefix, \
                 CREAL(vector->data.TYPE[lcv]), \
                 CIMAG(vector->data.TYPE[lcv])); \
        full = (strlen(str) > maxLength-2); \
        prefix = ","; \
    } \
    break;

    int lcv;
    char* prefix = "[";
    switch(vector->type.type) {
        APPEND_ELEMENTS_CASE(S8,char,"%hd")
        APPEND_ELEMENTS_CASE(S16,short int,"%hd")
        APPEND_ELEMENTS_CASE(S32,int,"%d")
        APPEND_ELEMENTS_CASE(S64,long,"%ld")
        APPEND_ELEMENTS_CASE(U8,unsigned char,"%hu")
        APPEND_ELEMENTS_CASE(U16,unsigned short,"%hu")
        APPEND_ELEMENTS_CASE(U32,unsigned int, "%u")
        APPEND_ELEMENTS_CASE(U64,unsigned long,"%lu")
        APPEND_ELEMENTS_CASE(F32,double,"%g")
        APPEND_ELEMENTS_CASE(F64,double,"%g")
        APPEND_ELEMENTS_CASE_COMPLEX(C32,crealf,cimagf)
        APPEND_ELEMENTS_CASE_COMPLEX(C64,creal,cimag)
    default:
        snprintf(str,maxLength,"[...]");
        break;
    }

    if (full) {
        // couldn't all fit in given string length

        // remove elements until there is room for ",...]"
        while (strlen(str) > maxLength - 5) {
            char* lastComma = strrchr(str,',');
            if (lastComma == NULL) { // no comma, must be first number
                str[1] = '\0';
            } else {
                *lastComma = '\0';
            }
        }
        strncat(str,",...]",maxLength);
    } else {
        strncat(str,"]",maxLength);
    }

    psFree(tempStr);

    return str;
}

psF64 p_psVectorGetElementF64(psVector* vector,
                              int position)
{
    if (vector == NULL) {
        return NAN;
    }
    if (position < 0 || position >= vector->n) {
        return NAN;
    }

    switch (vector->type.type) {
    case PS_TYPE_U8:
        return vector->data.U8[position];
        break;
    case PS_TYPE_U16:
        return vector->data.U16[position];
        break;
    case PS_TYPE_U32:
        return vector->data.U32[position];
        break;
    case PS_TYPE_U64:
        return vector->data.U64[position];
        break;
    case PS_TYPE_S8:
        return vector->data.S8[position];
        break;
    case PS_TYPE_S16:
        return vector->data.S16[position];
        break;
    case PS_TYPE_S32:
        return vector->data.S32[position];
        break;
    case PS_TYPE_S64:
        return vector->data.S64[position];
        break;
    case PS_TYPE_F32:
        return vector->data.F32[position];
        break;
    case PS_TYPE_F64:
        return vector->data.F64[position];
    default:
        return NAN;
    }
}

bool p_psVectorPrint (int fd,
                      psVector *a,
                      char *name)
{
    char line[1024];

    sprintf (line, "vector: %s\n", name);
    write (fd, line, strlen(line));

    for (int i = 0; i < a[0].n; i++) {
        sprintf (line, "%f\n", p_psVectorGetElementF64(a, i));
        write (fd, line, strlen(line));
    }
    sprintf (line, "\n");
    write (fd, line, strlen(line));
    return (true);
}


bool psVectorInit(psVector *vector,
                  ...)
{
    va_list argp;
    psU8  vU8;
    psU16 vU16;
    psU32 vU32;
    psU64 vU64;
    psS8 vS8;
    psS16 vS16;
    psS32 vS32;
    psS64 vS64;
    psF32 vF32;
    psF64 vF64;
    psC32 vC32;
    psC64 vC64;

    unsigned int temp;
    int temp2;

    if (vector == NULL)
        return (false);

    va_start (argp, vector);

    switch (vector->type.type) {
    case PS_TYPE_U8:
        temp = va_arg (argp, psU32);
        if ( temp >= PS_MIN_U8 && temp <= PS_MAX_U8 )
            vU8 = temp;
        else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error:  U8 Value out of Range.\n");
            return false;
        }
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.U8[iy] = vU8;
        }
        return (true);
    case PS_TYPE_F32:
        vF32 = (psF32)va_arg (argp, psF64);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.F32[iy] = vF32;
        }
        return (true);
    case PS_TYPE_F64:
        vF64 = va_arg (argp, psF64);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.F64[iy] = vF64;
        }
        return (true);
    case PS_TYPE_U16:
        temp = va_arg (argp, psU32);
        if ( temp >= PS_MIN_U16 && temp <= PS_MAX_U16 )
            vU16 = temp;
        else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error:  U16 Value out of Range.\n");
            return false;
        }
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.U16[iy] = vU16;
        }
        return (true);
    case PS_TYPE_U32:
        vU32 = va_arg (argp, psU32);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.U32[iy] = vU32;
        }
        return (true);
    case PS_TYPE_U64:
        vU64 = va_arg (argp, psU64);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.U64[iy] = vU64;
        }
        return (true);
    case PS_TYPE_S8:
        temp2 = va_arg (argp, psS32);
        if ( temp2 >= PS_MIN_S8 && temp2 <= PS_MAX_S8 )
            vS8 = temp2;
        else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error:  S8 Value out of Range.\n");
            return false;
        }
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.S8[iy] = vS8;
        }
        return (true);
    case PS_TYPE_S16:
        temp2 = va_arg (argp, psS32);
        if ( temp2 >= PS_MIN_S16 && temp2 <= PS_MAX_S16 )
            vS16 = temp2;
        else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error:  S16 Value out of Range.\n");
            return false;
        }
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.S16[iy] = vS16;
        }
        return (true);
    case PS_TYPE_S32:
        vS32 = va_arg (argp, psS32);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.S32[iy] = vS32;
        }
        return (true);
    case PS_TYPE_S64:
        vS64 = va_arg (argp, psS64);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.S64[iy] = vS64;
        }
        return (true);
    case PS_TYPE_C32:
        vC32 = va_arg (argp, psC32);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.C32[iy] = vC32;
        }
        return (true);
    case PS_TYPE_C64:
        vC64 = va_arg (argp, psC64);
        for (long iy = 0; iy < vector->n; iy++) {
            vector->data.C64[iy] = vC64;
        }
        return (true);
    default:
        psError (PS_ERR_BAD_PARAMETER_TYPE, true, "datatype %d not defined in psImageInit\n", vector->type);
        return (false);
    }
    return (false);
}

#define FUNC_MACRO_VECTOR_CREATE(TYPE) \
static psVector *vectorCreate##TYPE(psVector *input, double lower, double upper, double delta) \
{ \
    \
    int nBin = (upper - lower) / delta; \
    psVector *vec = psVectorRecycle(input, nBin, PS_TYPE_##TYPE); \
    for (int i = 0; i < nBin; i++) \
    { \
        vec->data.TYPE[i] = lower + (i * delta); \
    } \
    return vec; \
} \

FUNC_MACRO_VECTOR_CREATE(S8)
FUNC_MACRO_VECTOR_CREATE(S16)
FUNC_MACRO_VECTOR_CREATE(S32)
FUNC_MACRO_VECTOR_CREATE(S64)
FUNC_MACRO_VECTOR_CREATE(U8)
FUNC_MACRO_VECTOR_CREATE(U16)
FUNC_MACRO_VECTOR_CREATE(U32)
FUNC_MACRO_VECTOR_CREATE(U64)
FUNC_MACRO_VECTOR_CREATE(F32)
FUNC_MACRO_VECTOR_CREATE(F64)

psVector *psVectorCreate (psVector *input,
                          double lower,
                          double upper,
                          double delta,
                          psElemType type)
{
    psVector *out;
    switch (type) {
    case PS_TYPE_S8:
        out = vectorCreateS8(input, lower,  upper, delta);
        break;
    case PS_TYPE_S16:
        out = vectorCreateS16(input, lower, upper, delta);
        break;
    case PS_TYPE_S32:
        out = vectorCreateS32(input, lower, upper, delta);
        break;
    case PS_TYPE_S64:
        out = vectorCreateS64(input, lower, upper, delta);
        break;
    case PS_TYPE_U8:
        out = vectorCreateU8(input, lower, upper, delta);
        break;
    case PS_TYPE_U16:
        out = vectorCreateU16(input, lower, upper, delta);
        break;
    case PS_TYPE_U32:
        out = vectorCreateU32(input, lower, upper, delta);
        break;
    case PS_TYPE_U64:
        out = vectorCreateU64(input, lower, upper, delta);
        break;
    case PS_TYPE_F32:
        out = vectorCreateF32(input, lower, upper, delta);
        break;
    case PS_TYPE_F64:
        out = vectorCreateF64(input, lower, upper, delta);
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Invalid psType for Vector Create\n");
    }
    return (out);
}

bool psVectorSet(const psVector *input,
                 long position,
                 complex value)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psVector_NULL);
        return false;
    }
    if (position > input->n) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid position.  Number too large\n");
        return false;
    }
    if (position < 0)
        position += input->n;
    if (position < 0) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid position.  Negative number too large\n");
        return false;
    }

    if (position == input->n) {
        if (position >= input->nalloc) {
            psError(PS_ERR_BAD_PARAMETER_NULL, true,
                    "Specified position, %ld, is greater than n+1 of the vector, %ld.",
                    position, input->nalloc);
            return false;
        }
        (*(psVector**)&input)->n++;
    }

    switch (input->type.type) {
    case PS_TYPE_U8:
        *(psU8*)&input->data.U8[position] = value;
        break;
    case PS_TYPE_U16:
        *(psU16*)&input->data.U16[position] = value;
        break;
    case PS_TYPE_U32:
        *(psU32*)&input->data.U32[position] = value;
        break;
    case PS_TYPE_U64:
        *(psU64*)&input->data.U64[position] = value;
        break;
    case PS_TYPE_S8:
        *(psS8*)&input->data.S8[position] = value;
        break;
    case PS_TYPE_S16:
        *(psS16*)&input->data.S16[position] = value;
        break;
    case PS_TYPE_S32:
        *(psS32*)&input->data.S32[position] = value;
        break;
    case PS_TYPE_S64:
        *(psS64*)&input->data.S64[position] = value;
        break;
    case PS_TYPE_F32:
        *(psF32*)&input->data.F32[position] = value;
        break;
    case PS_TYPE_F64:
        *(psF64*)&input->data.F64[position] = value;
        break;
    case PS_TYPE_C32:
        *(psC32*)&input->data.C32[position] = value;
        break;
    case PS_TYPE_C64:
        *(psC64*)&input->data.C64[position] = value;
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Invalid psVector Data Type\n");
        return false;
    }

    return true;
}

complex psVectorGet(const psVector *input,
                    long position)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psVector_NULL);
        return NAN;
    }
    if (position >= input->n) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid position.  Number too large\n");
        return NAN;
    }
    if(position < 0)
        position += input->n;
    if (position < 0) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid position.  Negative number too large\n");
        return NAN;
    }
    switch (input->type.type) {
    case PS_TYPE_U8:
        return input->data.U8[position];
        break;
    case PS_TYPE_U16:
        return input->data.U16[position];
        break;
    case PS_TYPE_U32:
        return input->data.U32[position];
        break;
    case PS_TYPE_U64:
        return input->data.U64[position];
        break;
    case PS_TYPE_S8:
        return input->data.S8[position];
        break;
    case PS_TYPE_S16:
        return input->data.S16[position];
        break;
    case PS_TYPE_S32:
        return input->data.S32[position];
        break;
    case PS_TYPE_S64:
        return input->data.S64[position];
        break;
    case PS_TYPE_F32:
        return input->data.F32[position];
        break;
    case PS_TYPE_F64:
        return input->data.F64[position];
        break;
    case PS_TYPE_C32:
        return input->data.C32[position];
        break;
    case PS_TYPE_C64:
        return input->data.C64[position];
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Invalid psVector Data Type\n");
        return NAN;
    }

}

// count number of pixels with given mask value
long psVectorCountPixelMask (psVector *mask,
                             psMaskType value)
{
    long Npixels = 0;
    if (mask == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psVector_NULL);
        Npixels = -1;
        return Npixels;
    }

    psElemType type;
    type = mask->type.type;

    if (type != PS_TYPE_MASK) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                "psVector type does not match the specified psMaskType!\n");
        return -1;
    }

    switch (type) {
    case PS_TYPE_U8:
    case PS_TYPE_U16:
        for (long i = 0; i < mask->n; i++) {
            if (mask->data.PS_TYPE_MASK_DATA[i] & value) {
                Npixels ++;
            }
        }
        break;
    case PS_TYPE_S8:
    case PS_TYPE_S16:
    case PS_TYPE_S32:
    case PS_TYPE_S64:
    case PS_TYPE_U32:
    case PS_TYPE_U64:
    case PS_TYPE_F32:
    case PS_TYPE_F64:
    case PS_TYPE_C32:
    case PS_TYPE_C64:
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psVector_UNSUPPORTED_TYPE, type);
        return -1;
    }
    return (Npixels);
}

