IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1963


Ignore:
Timestamp:
Oct 5, 2004, 12:47:21 PM (22 years ago)
Author:
gusciora
Message:

Added macros in psConstants.h to check for NULL and mis-typed vectors.

Location:
trunk/psLib/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psConstants.h

    r1953 r1963  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-10-05 01:43:58 $
     8 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-10-05 22:47:21 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818#define RIGHT_SPLINE_DERIV 0.0
    1919
     20#define CONVERT_VECTOR_F64_TO_F32(OLD, NEW_PTR32, NEW_STATIC32) \
     21if (OLD->type.type == PS_TYPE_F32) { \
     22    NEW_PTR32 = (psVector *) OLD; \
     23} else if (OLD->type.type == PS_TYPE_F64) { \
     24    NEW_STATIC32 = psVectorRecycle(NEW_STATIC32, OLD->n, PS_TYPE_F32); \
     25    p_psMemSetPersistent(NEW_STATIC32, true); \
     26    p_psMemSetPersistent(NEW_STATIC32->data.V, true); \
     27    for (i=0; i < OLD->n ; i++) { \
     28        NEW_STATIC32->data.F32[i] = (float) OLD->data.F64[i]; \
     29    } \
     30    NEW_PTR32 = NEW_STATIC32; \
     31} \
     32
     33#define CONVERT_VECTOR_F32_TO_F64(OLD, NEW_PTR64, NEW_STATIC64) \
     34if (OLD->type.type == PS_TYPE_F64) { \
     35    NEW_PTR64 = (psVector *) OLD; \
     36} else if (OLD->type.type == PS_TYPE_F32) { \
     37    NEW_STATIC64 = psVectorRecycle(NEW_STATIC64, OLD->n, PS_TYPE_F64); \
     38    p_psMemSetPersistent(NEW_STATIC64, true); \
     39    p_psMemSetPersistent(NEW_STATIC64->data.V, true); \
     40    for (i=0; i < OLD->n ; i++) { \
     41        NEW_STATIC64->data.F64[i] = (double) OLD->data.F32[i]; \
     42    } \
     43    NEW_PTR64 = NEW_STATIC64; \
     44} \
     45
     46#define VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(VEC) \
     47if ((VEC->type.type != PS_TYPE_F32) && (VEC->type.type != PS_TYPE_F64)) { \
     48    psAbort(__func__, "Bad type for VEC (%d)", VEC->type.type); \
     49    return(NULL); \
     50} \
     51
     52#define GEN_PERSIST_YERR_VEC_F32(VEC, N) \
     53VEC = psVectorRecycle(VEC, N, PS_TYPE_F32); \
     54p_psMemSetPersistent(VEC, true); \
     55p_psMemSetPersistent(VEC->data.V, true); \
     56for (int i=0;i<N;i++) { \
     57    VEC->data.F32[i] = 1.0; \
     58} \
     59
     60#define GEN_PERSIST_YERR_VEC_F64(VEC, N) \
     61VEC = psVectorRecycle(VEC, N, PS_TYPE_F64); \
     62p_psMemSetPersistent(VEC, true); \
     63p_psMemSetPersistent(VEC->data.V, true); \
     64for (int i=0;i<N;i++) { \
     65    VEC->data.F64[i] = 1.0; \
     66} \
     67
     68#define GEN_PERSIST_X_INDEX_VEC_F32(VEC, N) \
     69VEC = psVectorRecycle(VEC, N, PS_TYPE_F32); \
     70p_psMemSetPersistent(VEC, true); \
     71p_psMemSetPersistent(VEC->data.V, true); \
     72for (int i=0;i<N;i++) { \
     73    VEC->data.F32[i] = (float) i; \
     74} \
     75
     76#define GEN_PERSIST_X_INDEX_VEC_F64(VEC, N) \
     77VEC = psVectorRecycle(VEC, N, PS_TYPE_F64); \
     78p_psMemSetPersistent(VEC, true); \
     79p_psMemSetPersistent(VEC->data.V, true); \
     80for (int i=0;i<N;i++) { \
     81    VEC->data.F64[i] = (float) i; \
     82} \
     83
    2084#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
    2185static psVector *NAME = NULL; \
    … …  
    2387p_psMemSetPersistent(NAME, true); \
    2488p_psMemSetPersistent(NAME->data.V, true); \
    25 
    26 
    27 
    2889
    2990/** Preprocessor macro to generate error on an incorrect type */
    … …  
    51112}
    52113
     114/** Preprocessor macro to generate error on a NULL vector */
     115#define PS_CHECK_NULL_VECTOR_RETURN_NULL(NAME) \
     116if (NAME == NULL || NAME->data.V == NULL) { \
     117    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
     118    return(NULL); \
     119} \
     120
    53121/** Preprocessor macro to generate error on a NULL poniter */
    54122#define PS_CHECK_NULL_PTR(NAME) \
    55123if (NAME == NULL) { \
    56124    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
     125}
     126
     127/** Preprocessor macro to generate error on a NULL poniter */
     128#define PS_CHECK_NULL_PTR_RETURN_NULL(NAME) \
     129if (NAME == NULL) { \
     130    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
     131    return(NULL); \
    57132}
    58133
    … …  
    75150}
    76151
     152/** Preprocessor macro to generate error for zero length vector */
     153#define PS_CHECK_EMPTY_VECTOR_RETURN_NULL(NAME) \
     154if (NAME->n < 1) { \
     155    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
     156    return(NULL); \
     157} \
     158
    77159/** Preprocessor macro to generate error on differing size vectors */
    78160#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \
    79161if (VEC1->n != VEC2->n) { \
    80162    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     163}
     164
     165/** Preprocessor macro to generate error on differing size vectors */
     166#define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(VEC1, VEC2) \
     167if (VEC1->n != VEC2->n) { \
     168    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     169    return(NULL); \
    81170}
    82171
    … …  
    112201}
    113202
     203#define PS_CHECK_NULL_1DPOLY_RETURN_NULL(NAME)                                                          \
     204if (NAME == NULL || NAME->coeff == NULL) {                                                         \
     205    psError(__func__,"Unallowable operation: %s or its coeffs is NULL.", #NAME);                          \
     206    return(NULL); \
     207} \
     208
    114209#define PS_PRINT_VECTOR(NAME) \
    115210for (int my_i=0;my_i<NAME->n;my_i++) { \

  • trunk/psLib/src/dataManip/psMinimize.c

    r1953 r1963  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-05 01:43:58 $
     11 *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-05 22:47:21 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    5656 
    5757XXX: change name
     58XXX: Use a static vector.
    5859 *****************************************************************************/
    5960void p_psBuildSums1D(double x,
    … …  
    6364    int i = 0;
    6465    double xSum = 0.0;
     66
    6567    if (sums == NULL) {
    6668        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
     69    }
     70
     71    if (polyOrder > sums->n) {
     72        sums = psVectorRealloc(sums, polyOrder);
    6773    }
    6874
    … …  
    7783VectorNormalizeGen(): this routine returns a psVector with "x" elements.  The
    7884values of the vector will be scaled uniformly between -1.0 and 1.0.
     85 
     86XXX: use a static vector.
    7987*****************************************************************************/
    8088psVector* VectorNormalizeGen(int x)
    … …  
    100108 
    101109XXX: This algorithm is derived from the Numerical Recipes.
     110XXX: use recyvled vectors for internal data.
    102111 *****************************************************************************/
    103112float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    … …  
    119128    float qn;
    120129
    121     if (x == NULL) {
    122         X = (float *) psAlloc(n * sizeof(float));
    123         for(i=0;i<n;i++) {
    124             X[i] = (float) i;
    125         }
    126         mustFreeX = true;
    127     }
    128 
    129130    // XXX: The second derivatives at the endpoints, undefined in the SDR,
    130131    // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV.
    … …  
    167168    return(derivs2);
    168169}
    169 
    170170
    171171/******************************************************************************
    … …  
    212212}
    213213
    214 
    215214/*****************************************************************************/
    216215/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    … …  
    258257    int i;
    259258    float slope;
    260 
    261     if (mySpline == NULL) {
    262         //XXX psErrorMsg()
    263     }
    264     if (y == NULL) {
    265         //XXX psErrorMsg()
    266     }
     259    psVector *x32 = NULL;
     260    psVector *y32 = NULL;
     261    psVector *yErr32 = NULL;
     262    static psVector *x32Static = NULL;
     263    static psVector *y32Static = NULL;
     264    static psVector *yErr32Static = NULL;
     265
     266    PS_CHECK_NULL_PTR_RETURN_NULL(mySpline);
     267    PS_CHECK_NULL_PTR_RETURN_NULL(y);
     268    VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(y);
     269    CONVERT_VECTOR_F64_TO_F32(y, y32, y32Static);
     270
     271    // If yErr==NULL, set all errors equal.
     272    if (yErr == NULL) {
     273        GEN_PERSIST_YERR_VEC_F32(yErr32Static, y->n);
     274        yErr32 = yErr32Static;
     275    } else {
     276        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(yErr);
     277        CONVERT_VECTOR_F64_TO_F32(yErr, yErr32, yErr32Static);
     278    }
     279
     280    // If x==NULL, create an x32 vector with x values set to (0:n).
     281    if (x == NULL) {
     282        GEN_PERSIST_X_INDEX_VEC_F32(x32Static, y->n);
     283        x32 = x32Static;
     284    } else {
     285        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(x);
     286        CONVERT_VECTOR_F64_TO_F32(x, x32, x32Static);
     287    }
     288    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y);
     289    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr, y);
    267290
    268291    if (y->n != (1 + mySpline->n)) {
    269         psLogMsg(__func__, PS_LOG_WARN,
    270                  "data size / spline size mismatch (%d %d)\n",
    271                  y->n, mySpline->n);
    272         // XXX: psErrorMsg()
     292        psError(__func__, "data size / spline size mismatch (%d %d)\n",
     293                y->n, mySpline->n);
     294        return(NULL);
    273295    }
    274296
    … …  
    295317    // Check if these are cubic splines (n==4).  If not, psError.
    296318    if (4 != (mySpline->spline[0])->n) {
    297         psLogMsg(__func__, PS_LOG_WARN,
    298                  "Don't know how to generate %d-order splines.",
    299                  (mySpline->spline[0])->n-1);
    300         // XXX: psErrorMsg()
     319        psError(__func__, "Don't know how to generate %d-order splines.",
     320                (mySpline->spline[0])->n-1);
    301321        return(NULL);
    302322    }
    … …  
    315335        H = x->data.F32[i+1] - x->data.F32[i];
    316336        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    317                 "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
     337                "x data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
    318338        //
    319339        // ******** Calculate 0-order term ********
    … …  
    630650}
    631651
    632 #define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
    633 static psVector *NAME = NULL; \
    634 NAME = psVectorRecycle(NAME, SIZE, TYPE); \
    635 p_psMemSetPersistent(NAME, true); \
    636 p_psMemSetPersistent(NAME->data.V, true); \
    637 
    638 
    639652/******************************************************************************
    640653p_psVectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
    … …  
    854867    int i;
    855868    psPolynomial1D *tmpPoly;
    856     PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
    857     PS_CHECK_NULL_VECTOR_ACTION(y, 0);
    858     PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
     869    PS_CHECK_NULL_1DPOLY_RETURN_NULL(myPoly);
     870    PS_CHECK_NULL_VECTOR_RETURN_NULL(y);
     871    PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y);
    859872    psVector *x64 = NULL;
    860873    psVector *y64 = NULL;
    … …  
    864877    static psVector *yErr64Static = NULL;
    865878
    866     if (y->type.type == PS_TYPE_F64) {
    867         y64 = (psVector *) y;
    868     } else if (y->type.type == PS_TYPE_F32) {
    869         y64Static = psVectorRecycle(y64Static, y->n, PS_TYPE_F64);
    870         p_psMemSetPersistent(y64Static, true);
    871         p_psMemSetPersistent(y64Static->data.V, true);
    872         y64 = y64Static;
    873     } else {
    874         // XXX: psError() bad type.
    875         psAbort(__func__, "Bad type for y64 (%d)", y->type.type);
    876     }
    877 
     879    VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(y);
     880    CONVERT_VECTOR_F32_TO_F64(y, y64, y64Static);
    878881
    879882    // If yErr==NULL, set all errors equal.
    880883    if (yErr == NULL) {
    881         yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64);
    882         p_psMemSetPersistent(yErr64Static, true);
    883         p_psMemSetPersistent(yErr64Static->data.V, true);
    884 
    885         for (i=0;i<yErr64Static->n;i++) {
    886             yErr64Static->data.F64[i] = 1.0;
    887         }
     884        GEN_PERSIST_YERR_VEC_F64(yErr64Static, y->n);
    888885        yErr64 = yErr64Static;
    889886    } else {
    890         if (yErr->type.type == PS_TYPE_F64) {
    891             yErr64 = (psVector *) yErr;
    892         } else if (yErr->type.type == PS_TYPE_F32) {
    893             yErr64Static = psVectorRecycle(yErr64Static, yErr->n, PS_TYPE_F64);
    894             p_psMemSetPersistent(yErr64Static, true);
    895             p_psMemSetPersistent(yErr64Static->data.V, true);
    896         } else {
    897             // XXX: psError() bad type.
    898             psAbort(__func__, "Bad type for yErr64");
    899         }
    900     }
    901 
    902     // If x==NULL, create an x64 vector with x values set to (0:n), and if
    903     // this is a Chebyshev polynomial, we must scale to (-1:1).
    904     // XXX: Verify that this is the correct action.
     887        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(yErr);
     888        CONVERT_VECTOR_F32_TO_F64(yErr, yErr64, yErr64Static);
     889    }
     890
     891    // If x==NULL, create an x64 vector with x values set to (0:n).
    905892    if (x == NULL) {
    906         x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64);
    907         p_psMemSetPersistent(x64Static, true);
    908         p_psMemSetPersistent(x64Static->data.V, true);
     893        GEN_PERSIST_X_INDEX_VEC_F64(x64Static, y->n);
     894        if (myPoly->type == PS_POLYNOMIAL_CHEB) {
     895            p_psNormalizeVector(x64Static);
     896        }
    909897        x64 = x64Static;
    910 
    911         if (myPoly->type == PS_POLYNOMIAL_ORD) {
    912             for (i=0;i<x64->n;i++) {
    913                 x64->data.F64[i] = (float) i;
    914             }
    915         } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    916             double min = 0.0;
    917             double max = (double) (y->n - 1);
    918 
    919             for (i=0;i<x64->n;i++) {
    920                 x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) /
    921                                    (0.5 * (max - min));
    922             }
    923         }
    924898    } else {
    925         if (x->type.type == PS_TYPE_F64) {
    926             x64 = (psVector *) x;
    927         } else if (x->type.type == PS_TYPE_F32) {
    928             x64Static = psVectorRecycle(x64Static, x->n, PS_TYPE_F64);
    929             p_psMemSetPersistent(x64Static, true);
    930             p_psMemSetPersistent(x64Static->data.V, true);
    931             x64 = x64Static;
    932         } else {
    933             // XXX: psError() bad type.
    934             psAbort(__func__, "Bad type for x64");
    935         }
    936     }
    937 
    938     PS_CHECK_VECTOR_SIZE_EQUAL(y64, x64);
    939     PS_CHECK_VECTOR_SIZE_EQUAL(y64, yErr64);
     899        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(x);
     900        CONVERT_VECTOR_F32_TO_F64(x, x64, x64Static);
     901    }
     902    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x64, y64);
     903    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr64, y64);
    940904
    941905    // Call the appropriate vector fitting routine.
    942906    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    943         tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     907        p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    944908    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    945909        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    … …  
    948912        return(NULL);
    949913    }
    950 
    951     // Free any allocated memory.
    952     /*
    953         if ((yErr == NULL) || (yErr->type.type != PS_TYPE_F64)) {
    954             psFree(yErr64);
    955         }
    956         if ((x == NULL) || (x->type.type != PS_TYPE_F64)) {
    957             psFree(x64);
    958         }
    959         if ((y == NULL) || (y->type.type != PS_TYPE_F64)) {
    960             psFree(y64);
    961         }
    962     */
    963914
    964915    return(myPoly);

  • trunk/psLib/src/dataManip/psMinimize.h

    r1953 r1963  
    88 *  @author George Gusciora, MHPCC
    99 *
    10  *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-10-05 01:43:58 $
     10 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-10-05 22:47:21 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2828#include "psMatrix.h"
    2929#include "psFunctions.h"
     30#include "psStats.h"
    3031#include "psTrace.h"
    3132#include "psLogMsg.h"

  • trunk/psLib/src/dataManip/psStats.c

    r1928 r1963  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-29 19:52:53 $
     11 *  @version $Revision: 1.64 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-05 22:47:21 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    836836
    837837/*****************************************************************************
    838 p_psNormalizeVector(myData): this is a private function which normalizes the
     838p_psNormalizeVectorF32(myData): this is a private function which normalizes the
    839839elements of a vector to a range between 0.0 and 1.0.
    840840 *****************************************************************************/
    841 void p_psNormalizeVector(psVector* myData)
     841void p_psNormalizeVectorF32(psVector* myData)
    842842{
    843843    float min = (float)HUGE;
    … …  
    867867p_psNormalizeVectorF64(myData): this is a private function which normalizes the
    868868elements of a vector to a range between -1.0 and 1.0.
    869  
    870 XXX: incorporate this into above
    871  
    872869XXX: 0-1 or -1:1?
    873870 *****************************************************************************/
    874871void p_psNormalizeVectorF64(psVector* myData)
    875872{
    876     float min = (float)HUGE;
    877     float max = (float)-HUGE;
    878     float range = 0.0;
     873    float min = (double)HUGE;
     874    float max = (double)-HUGE;
     875    double range = 0.0;
    879876    int i = 0;
    880877
    … …  
    892889        myData->data.F64[i] = -1.0 + 2.0 *
    893890                              ((myData->data.F64[i] - min) / range);
     891    }
     892
     893    //    for (i = 0; i < myData->n; i++) {
     894    //        myData->data.F64[i] = (myData->data.F64[i] - 0.5 * (min + max)) /
     895    //                              (0.5 * (max - min));
     896    //    }
     897}
     898
     899void p_psNormalizeVector(psVector* myData)
     900{
     901    if (myData->type.type == PS_TYPE_F32) {
     902        p_psNormalizeVectorF32(myData);
     903    } else if (myData->type.type == PS_TYPE_F64) {
     904        p_psNormalizeVectorF64(myData);
     905    } else {
     906        psError(__func__, "Unalowable data type.\n");
    894907    }
    895908}
    … …  
    11911204    // The following was necessary to fit a gaussian to the data, since
    11921205    // gaussian functions produce data between 0.0 and 1.0.
    1193     // p_psNormalizeVector(robustHistogramVector);
     1206    // p_psNormalizeVectorF32(robustHistogramVector);
    11941207
    11951208    /**************************************************************************

  • trunk/psLib/src/dataManip/psStats.h

    r1907 r1963  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-09-27 23:41:42 $
     12 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-10-05 22:47:21 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    176176
    177177void p_psNormalizeVector(psVector* myData);
    178 void p_psNormalizeVectorF64(psVector* myData);
    179 
    180178
    181179/// @}

  • trunk/psLib/src/math/psConstants.h

    r1953 r1963  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-10-05 01:43:58 $
     8 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-10-05 22:47:21 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818#define RIGHT_SPLINE_DERIV 0.0
    1919
     20#define CONVERT_VECTOR_F64_TO_F32(OLD, NEW_PTR32, NEW_STATIC32) \
     21if (OLD->type.type == PS_TYPE_F32) { \
     22    NEW_PTR32 = (psVector *) OLD; \
     23} else if (OLD->type.type == PS_TYPE_F64) { \
     24    NEW_STATIC32 = psVectorRecycle(NEW_STATIC32, OLD->n, PS_TYPE_F32); \
     25    p_psMemSetPersistent(NEW_STATIC32, true); \
     26    p_psMemSetPersistent(NEW_STATIC32->data.V, true); \
     27    for (i=0; i < OLD->n ; i++) { \
     28        NEW_STATIC32->data.F32[i] = (float) OLD->data.F64[i]; \
     29    } \
     30    NEW_PTR32 = NEW_STATIC32; \
     31} \
     32
     33#define CONVERT_VECTOR_F32_TO_F64(OLD, NEW_PTR64, NEW_STATIC64) \
     34if (OLD->type.type == PS_TYPE_F64) { \
     35    NEW_PTR64 = (psVector *) OLD; \
     36} else if (OLD->type.type == PS_TYPE_F32) { \
     37    NEW_STATIC64 = psVectorRecycle(NEW_STATIC64, OLD->n, PS_TYPE_F64); \
     38    p_psMemSetPersistent(NEW_STATIC64, true); \
     39    p_psMemSetPersistent(NEW_STATIC64->data.V, true); \
     40    for (i=0; i < OLD->n ; i++) { \
     41        NEW_STATIC64->data.F64[i] = (double) OLD->data.F32[i]; \
     42    } \
     43    NEW_PTR64 = NEW_STATIC64; \
     44} \
     45
     46#define VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(VEC) \
     47if ((VEC->type.type != PS_TYPE_F32) && (VEC->type.type != PS_TYPE_F64)) { \
     48    psAbort(__func__, "Bad type for VEC (%d)", VEC->type.type); \
     49    return(NULL); \
     50} \
     51
     52#define GEN_PERSIST_YERR_VEC_F32(VEC, N) \
     53VEC = psVectorRecycle(VEC, N, PS_TYPE_F32); \
     54p_psMemSetPersistent(VEC, true); \
     55p_psMemSetPersistent(VEC->data.V, true); \
     56for (int i=0;i<N;i++) { \
     57    VEC->data.F32[i] = 1.0; \
     58} \
     59
     60#define GEN_PERSIST_YERR_VEC_F64(VEC, N) \
     61VEC = psVectorRecycle(VEC, N, PS_TYPE_F64); \
     62p_psMemSetPersistent(VEC, true); \
     63p_psMemSetPersistent(VEC->data.V, true); \
     64for (int i=0;i<N;i++) { \
     65    VEC->data.F64[i] = 1.0; \
     66} \
     67
     68#define GEN_PERSIST_X_INDEX_VEC_F32(VEC, N) \
     69VEC = psVectorRecycle(VEC, N, PS_TYPE_F32); \
     70p_psMemSetPersistent(VEC, true); \
     71p_psMemSetPersistent(VEC->data.V, true); \
     72for (int i=0;i<N;i++) { \
     73    VEC->data.F32[i] = (float) i; \
     74} \
     75
     76#define GEN_PERSIST_X_INDEX_VEC_F64(VEC, N) \
     77VEC = psVectorRecycle(VEC, N, PS_TYPE_F64); \
     78p_psMemSetPersistent(VEC, true); \
     79p_psMemSetPersistent(VEC->data.V, true); \
     80for (int i=0;i<N;i++) { \
     81    VEC->data.F64[i] = (float) i; \
     82} \
     83
    2084#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
    2185static psVector *NAME = NULL; \
    … …  
    2387p_psMemSetPersistent(NAME, true); \
    2488p_psMemSetPersistent(NAME->data.V, true); \
    25 
    26 
    27 
    2889
    2990/** Preprocessor macro to generate error on an incorrect type */
    … …  
    51112}
    52113
     114/** Preprocessor macro to generate error on a NULL vector */
     115#define PS_CHECK_NULL_VECTOR_RETURN_NULL(NAME) \
     116if (NAME == NULL || NAME->data.V == NULL) { \
     117    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
     118    return(NULL); \
     119} \
     120
    53121/** Preprocessor macro to generate error on a NULL poniter */
    54122#define PS_CHECK_NULL_PTR(NAME) \
    55123if (NAME == NULL) { \
    56124    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
     125}
     126
     127/** Preprocessor macro to generate error on a NULL poniter */
     128#define PS_CHECK_NULL_PTR_RETURN_NULL(NAME) \
     129if (NAME == NULL) { \
     130    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
     131    return(NULL); \
    57132}
    58133
    … …  
    75150}
    76151
     152/** Preprocessor macro to generate error for zero length vector */
     153#define PS_CHECK_EMPTY_VECTOR_RETURN_NULL(NAME) \
     154if (NAME->n < 1) { \
     155    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
     156    return(NULL); \
     157} \
     158
    77159/** Preprocessor macro to generate error on differing size vectors */
    78160#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \
    79161if (VEC1->n != VEC2->n) { \
    80162    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     163}
     164
     165/** Preprocessor macro to generate error on differing size vectors */
     166#define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(VEC1, VEC2) \
     167if (VEC1->n != VEC2->n) { \
     168    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     169    return(NULL); \
    81170}
    82171
    … …  
    112201}
    113202
     203#define PS_CHECK_NULL_1DPOLY_RETURN_NULL(NAME)                                                          \
     204if (NAME == NULL || NAME->coeff == NULL) {                                                         \
     205    psError(__func__,"Unallowable operation: %s or its coeffs is NULL.", #NAME);                          \
     206    return(NULL); \
     207} \
     208
    114209#define PS_PRINT_VECTOR(NAME) \
    115210for (int my_i=0;my_i<NAME->n;my_i++) { \

  • trunk/psLib/src/math/psMinimize.c

    r1953 r1963  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-05 01:43:58 $
     11 *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-05 22:47:21 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    5656 
    5757XXX: change name
     58XXX: Use a static vector.
    5859 *****************************************************************************/
    5960void p_psBuildSums1D(double x,
    … …  
    6364    int i = 0;
    6465    double xSum = 0.0;
     66
    6567    if (sums == NULL) {
    6668        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
     69    }
     70
     71    if (polyOrder > sums->n) {
     72        sums = psVectorRealloc(sums, polyOrder);
    6773    }
    6874
    … …  
    7783VectorNormalizeGen(): this routine returns a psVector with "x" elements.  The
    7884values of the vector will be scaled uniformly between -1.0 and 1.0.
     85 
     86XXX: use a static vector.
    7987*****************************************************************************/
    8088psVector* VectorNormalizeGen(int x)
    … …  
    100108 
    101109XXX: This algorithm is derived from the Numerical Recipes.
     110XXX: use recyvled vectors for internal data.
    102111 *****************************************************************************/
    103112float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    … …  
    119128    float qn;
    120129
    121     if (x == NULL) {
    122         X = (float *) psAlloc(n * sizeof(float));
    123         for(i=0;i<n;i++) {
    124             X[i] = (float) i;
    125         }
    126         mustFreeX = true;
    127     }
    128 
    129130    // XXX: The second derivatives at the endpoints, undefined in the SDR,
    130131    // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV.
    … …  
    167168    return(derivs2);
    168169}
    169 
    170170
    171171/******************************************************************************
    … …  
    212212}
    213213
    214 
    215214/*****************************************************************************/
    216215/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    … …  
    258257    int i;
    259258    float slope;
    260 
    261     if (mySpline == NULL) {
    262         //XXX psErrorMsg()
    263     }
    264     if (y == NULL) {
    265         //XXX psErrorMsg()
    266     }
     259    psVector *x32 = NULL;
     260    psVector *y32 = NULL;
     261    psVector *yErr32 = NULL;
     262    static psVector *x32Static = NULL;
     263    static psVector *y32Static = NULL;
     264    static psVector *yErr32Static = NULL;
     265
     266    PS_CHECK_NULL_PTR_RETURN_NULL(mySpline);
     267    PS_CHECK_NULL_PTR_RETURN_NULL(y);
     268    VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(y);
     269    CONVERT_VECTOR_F64_TO_F32(y, y32, y32Static);
     270
     271    // If yErr==NULL, set all errors equal.
     272    if (yErr == NULL) {
     273        GEN_PERSIST_YERR_VEC_F32(yErr32Static, y->n);
     274        yErr32 = yErr32Static;
     275    } else {
     276        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(yErr);
     277        CONVERT_VECTOR_F64_TO_F32(yErr, yErr32, yErr32Static);
     278    }
     279
     280    // If x==NULL, create an x32 vector with x values set to (0:n).
     281    if (x == NULL) {
     282        GEN_PERSIST_X_INDEX_VEC_F32(x32Static, y->n);
     283        x32 = x32Static;
     284    } else {
     285        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(x);
     286        CONVERT_VECTOR_F64_TO_F32(x, x32, x32Static);
     287    }
     288    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y);
     289    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr, y);
    267290
    268291    if (y->n != (1 + mySpline->n)) {
    269         psLogMsg(__func__, PS_LOG_WARN,
    270                  "data size / spline size mismatch (%d %d)\n",
    271                  y->n, mySpline->n);
    272         // XXX: psErrorMsg()
     292        psError(__func__, "data size / spline size mismatch (%d %d)\n",
     293                y->n, mySpline->n);
     294        return(NULL);
    273295    }
    274296
    … …  
    295317    // Check if these are cubic splines (n==4).  If not, psError.
    296318    if (4 != (mySpline->spline[0])->n) {
    297         psLogMsg(__func__, PS_LOG_WARN,
    298                  "Don't know how to generate %d-order splines.",
    299                  (mySpline->spline[0])->n-1);
    300         // XXX: psErrorMsg()
     319        psError(__func__, "Don't know how to generate %d-order splines.",
     320                (mySpline->spline[0])->n-1);
    301321        return(NULL);
    302322    }
    … …  
    315335        H = x->data.F32[i+1] - x->data.F32[i];
    316336        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    317                 "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
     337                "x data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
    318338        //
    319339        // ******** Calculate 0-order term ********
    … …  
    630650}
    631651
    632 #define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
    633 static psVector *NAME = NULL; \
    634 NAME = psVectorRecycle(NAME, SIZE, TYPE); \
    635 p_psMemSetPersistent(NAME, true); \
    636 p_psMemSetPersistent(NAME->data.V, true); \
    637 
    638 
    639652/******************************************************************************
    640653p_psVectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
    … …  
    854867    int i;
    855868    psPolynomial1D *tmpPoly;
    856     PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
    857     PS_CHECK_NULL_VECTOR_ACTION(y, 0);
    858     PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
     869    PS_CHECK_NULL_1DPOLY_RETURN_NULL(myPoly);
     870    PS_CHECK_NULL_VECTOR_RETURN_NULL(y);
     871    PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y);
    859872    psVector *x64 = NULL;
    860873    psVector *y64 = NULL;
    … …  
    864877    static psVector *yErr64Static = NULL;
    865878
    866     if (y->type.type == PS_TYPE_F64) {
    867         y64 = (psVector *) y;
    868     } else if (y->type.type == PS_TYPE_F32) {
    869         y64Static = psVectorRecycle(y64Static, y->n, PS_TYPE_F64);
    870         p_psMemSetPersistent(y64Static, true);
    871         p_psMemSetPersistent(y64Static->data.V, true);
    872         y64 = y64Static;
    873     } else {
    874         // XXX: psError() bad type.
    875         psAbort(__func__, "Bad type for y64 (%d)", y->type.type);
    876     }
    877 
     879    VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(y);
     880    CONVERT_VECTOR_F32_TO_F64(y, y64, y64Static);
    878881
    879882    // If yErr==NULL, set all errors equal.
    880883    if (yErr == NULL) {
    881         yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64);
    882         p_psMemSetPersistent(yErr64Static, true);
    883         p_psMemSetPersistent(yErr64Static->data.V, true);
    884 
    885         for (i=0;i<yErr64Static->n;i++) {
    886             yErr64Static->data.F64[i] = 1.0;
    887         }
     884        GEN_PERSIST_YERR_VEC_F64(yErr64Static, y->n);
    888885        yErr64 = yErr64Static;
    889886    } else {
    890         if (yErr->type.type == PS_TYPE_F64) {
    891             yErr64 = (psVector *) yErr;
    892         } else if (yErr->type.type == PS_TYPE_F32) {
    893             yErr64Static = psVectorRecycle(yErr64Static, yErr->n, PS_TYPE_F64);
    894             p_psMemSetPersistent(yErr64Static, true);
    895             p_psMemSetPersistent(yErr64Static->data.V, true);
    896         } else {
    897             // XXX: psError() bad type.
    898             psAbort(__func__, "Bad type for yErr64");
    899         }
    900     }
    901 
    902     // If x==NULL, create an x64 vector with x values set to (0:n), and if
    903     // this is a Chebyshev polynomial, we must scale to (-1:1).
    904     // XXX: Verify that this is the correct action.
     887        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(yErr);
     888        CONVERT_VECTOR_F32_TO_F64(yErr, yErr64, yErr64Static);
     889    }
     890
     891    // If x==NULL, create an x64 vector with x values set to (0:n).
    905892    if (x == NULL) {
    906         x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64);
    907         p_psMemSetPersistent(x64Static, true);
    908         p_psMemSetPersistent(x64Static->data.V, true);
     893        GEN_PERSIST_X_INDEX_VEC_F64(x64Static, y->n);
     894        if (myPoly->type == PS_POLYNOMIAL_CHEB) {
     895            p_psNormalizeVector(x64Static);
     896        }
    909897        x64 = x64Static;
    910 
    911         if (myPoly->type == PS_POLYNOMIAL_ORD) {
    912             for (i=0;i<x64->n;i++) {
    913                 x64->data.F64[i] = (float) i;
    914             }
    915         } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    916             double min = 0.0;
    917             double max = (double) (y->n - 1);
    918 
    919             for (i=0;i<x64->n;i++) {
    920                 x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) /
    921                                    (0.5 * (max - min));
    922             }
    923         }
    924898    } else {
    925         if (x->type.type == PS_TYPE_F64) {
    926             x64 = (psVector *) x;
    927         } else if (x->type.type == PS_TYPE_F32) {
    928             x64Static = psVectorRecycle(x64Static, x->n, PS_TYPE_F64);
    929             p_psMemSetPersistent(x64Static, true);
    930             p_psMemSetPersistent(x64Static->data.V, true);
    931             x64 = x64Static;
    932         } else {
    933             // XXX: psError() bad type.
    934             psAbort(__func__, "Bad type for x64");
    935         }
    936     }
    937 
    938     PS_CHECK_VECTOR_SIZE_EQUAL(y64, x64);
    939     PS_CHECK_VECTOR_SIZE_EQUAL(y64, yErr64);
     899        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(x);
     900        CONVERT_VECTOR_F32_TO_F64(x, x64, x64Static);
     901    }
     902    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x64, y64);
     903    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr64, y64);
    940904
    941905    // Call the appropriate vector fitting routine.
    942906    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    943         tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     907        p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    944908    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    945909        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    … …  
    948912        return(NULL);
    949913    }
    950 
    951     // Free any allocated memory.
    952     /*
    953         if ((yErr == NULL) || (yErr->type.type != PS_TYPE_F64)) {
    954             psFree(yErr64);
    955         }
    956         if ((x == NULL) || (x->type.type != PS_TYPE_F64)) {
    957             psFree(x64);
    958         }
    959         if ((y == NULL) || (y->type.type != PS_TYPE_F64)) {
    960             psFree(y64);
    961         }
    962     */
    963914
    964915    return(myPoly);

  • trunk/psLib/src/math/psMinimize.h

    r1953 r1963  
    88 *  @author George Gusciora, MHPCC
    99 *
    10  *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-10-05 01:43:58 $
     10 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-10-05 22:47:21 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2828#include "psMatrix.h"
    2929#include "psFunctions.h"
     30#include "psStats.h"
    3031#include "psTrace.h"
    3132#include "psLogMsg.h"

  • trunk/psLib/src/math/psStats.c

    r1928 r1963  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-29 19:52:53 $
     11 *  @version $Revision: 1.64 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-05 22:47:21 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    836836
    837837/*****************************************************************************
    838 p_psNormalizeVector(myData): this is a private function which normalizes the
     838p_psNormalizeVectorF32(myData): this is a private function which normalizes the
    839839elements of a vector to a range between 0.0 and 1.0.
    840840 *****************************************************************************/
    841 void p_psNormalizeVector(psVector* myData)
     841void p_psNormalizeVectorF32(psVector* myData)
    842842{
    843843    float min = (float)HUGE;
    … …  
    867867p_psNormalizeVectorF64(myData): this is a private function which normalizes the
    868868elements of a vector to a range between -1.0 and 1.0.
    869  
    870 XXX: incorporate this into above
    871  
    872869XXX: 0-1 or -1:1?
    873870 *****************************************************************************/
    874871void p_psNormalizeVectorF64(psVector* myData)
    875872{
    876     float min = (float)HUGE;
    877     float max = (float)-HUGE;
    878     float range = 0.0;
     873    float min = (double)HUGE;
     874    float max = (double)-HUGE;
     875    double range = 0.0;
    879876    int i = 0;
    880877
    … …  
    892889        myData->data.F64[i] = -1.0 + 2.0 *
    893890                              ((myData->data.F64[i] - min) / range);
     891    }
     892
     893    //    for (i = 0; i < myData->n; i++) {
     894    //        myData->data.F64[i] = (myData->data.F64[i] - 0.5 * (min + max)) /
     895    //                              (0.5 * (max - min));
     896    //    }
     897}
     898
     899void p_psNormalizeVector(psVector* myData)
     900{
     901    if (myData->type.type == PS_TYPE_F32) {
     902        p_psNormalizeVectorF32(myData);
     903    } else if (myData->type.type == PS_TYPE_F64) {
     904        p_psNormalizeVectorF64(myData);
     905    } else {
     906        psError(__func__, "Unalowable data type.\n");
    894907    }
    895908}
    … …  
    11911204    // The following was necessary to fit a gaussian to the data, since
    11921205    // gaussian functions produce data between 0.0 and 1.0.
    1193     // p_psNormalizeVector(robustHistogramVector);
     1206    // p_psNormalizeVectorF32(robustHistogramVector);
    11941207
    11951208    /**************************************************************************

  • trunk/psLib/src/math/psStats.h

    r1907 r1963  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-09-27 23:41:42 $
     12 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-10-05 22:47:21 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    176176
    177177void p_psNormalizeVector(psVector* myData);
    178 void p_psNormalizeVectorF64(psVector* myData);
    179 
    180178
    181179/// @}
Note: See TracChangeset for help on using the changeset viewer.