IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2908


Ignore:
Timestamp:
Jan 5, 2005, 10:59:03 AM (22 years ago)
Author:
gusciora
Message:

* empty log message *

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/dataManip/psMinimize.c

    r2788 r2908  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-22 05:09:32 $
     11 *  @version $Revision: 1.97 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-01-05 20:59:03 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    190190    n = spline->n;
    191191    klo = p_psVectorBinDisect32(spline->domains, (spline->n)+1, X);
     192    if (klo < 0) {
     193        psLogMsg(__func__, PS_LOG_WARN,
     194                 "WARNING: psMinimize.c: p_psNRSpline1DEval(): p_psVectorBinDisect32 returned an error (%d).\n", klo);
     195        return(NAN);
     196    }
    192197    khi = klo + 1;
    193198    H = (spline->domains)[khi] - (spline->domains)[klo];
     
    245250                                const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
    246251{
    247     PS_PTR_CHECK_NULL(mySpline, NULL);
    248     PS_INT_CHECK_NON_NEGATIVE(mySpline->n, NULL);
    249     PS_VECTOR_CHECK_NULL(x, NULL);
    250252    PS_VECTOR_CHECK_NULL(y, NULL);
    251253    PS_VECTOR_CHECK_TYPE_F32_OR_F64(y, NULL);
     
    278280    // If x==NULL, create an x32 vector with x values set to (0:n).
    279281    if (x == NULL) {
     282
    280283        PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n);
    281284        x32 = x32Static;
     
    286289    PS_VECTOR_CHECK_SIZE_EQUAL(x32, y32, NULL);
    287290    PS_VECTOR_CHECK_SIZE_EQUAL(yErr32, y32, NULL);
     291
     292    /*
     293        XXX:
     294        This can not be implemented until SDR states what order spline should be
     295        created.
     296     
     297        Should we error if mySPline is not NULL?
     298    */
     299    if (mySpline == NULL) {
     300        mySpline = psSpline1DAllocGeneric(x32, 3);
     301    }
     302    PS_PTR_CHECK_NULL(mySpline, NULL);
     303    PS_INT_CHECK_NON_NEGATIVE(mySpline->n, NULL);
     304
    288305
    289306    if (y32->n != (1 + mySpline->n)) {
     
    335352        H = x32->data.F32[i+1] - x32->data.F32[i];
    336353        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    337                 "x data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
     354                "x data (%f - %f) (%f)\n",
     355                x32->data.F32[i],
     356                x32->data.F32[i+1], H);
    338357        //
    339358        // ******** Calculate 0-order term ********
     
    429448        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
    430449    } else {
    431         // XXX: Check size of derivative
     450        PS_IMAGE_CHECK_SIZE(deriv, params->n, coords->n, NULL);
     451        PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    432452    }
    433453
     
    499519    } else {
    500520        PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL);
     521        PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    501522    }
    502523
     
    715736        // Solve A * alpha = Beta
    716737        //
     738        // XXX: How do we know if these functions were successful?
     739        //
    717740        aOut = psMatrixLUD(aOut, perm, A);
    718741        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
     
    963986    }
    964987
     988    // XXX: How do we know if these routines were successful?
    965989    ALUD = psMatrixLUD(ALUD, outPerm, A);
    966990    coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
     
    970994        // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]);
    971995    }
    972 
    973     // for (i=0;i<x->n;i++) {
    974     // printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i],
    975     // myPoly));
    976     // }
    977996
    978997    psFree(A);
     
    10081027
    10091028    psS32 i;
    1010     psPolynomial1D *tmpPoly;
    10111029    psVector *x64 = NULL;
    10121030    psVector *y64 = NULL;
     
    10441062
    10451063    // Call the appropriate vector fitting routine.
     1064    psPolynomial1D *rc = NULL;
    10461065    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    1047         VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     1066        rc = VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    10481067    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    1049         tmpPoly = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
     1068        rc = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    10501069    } else {
    10511070        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    10521071                "unknown polynomial type.\n");
     1072        return(NULL);
     1073    }
     1074    if (rc == NULL) {
     1075        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
    10531076        return(NULL);
    10541077    }
     
    14481471    if (bracket == NULL) {
    14491472        psError(PS_ERR_UNKNOWN, false,
    1450                 "Could not bracket minimum.");
     1473                "Could not bracket minimum.  Returning NAN.\n");
    14511474        return(NAN);
    14521475    }
     
    16361659                if (isnan(mul)) {
    16371660                    psError(PS_ERR_UNKNOWN, false,
    1638                             "Could not perform line minimization");
     1661                            "Could not perform line minimization.  Returning FALSE.\n");
    16391662                    psFree(v);
    16401663                    return(false);
     
    16771700        if (isnan(mul)) {
    16781701            psError(PS_ERR_UNKNOWN, false,
    1679                     "Could not perform line minimization.");
     1702                    "Could not perform line minimization.  Returning FALSE.\n");
    16801703            psFree(v);
    16811704            return(false);
  • trunk/psLib/src/math/psMinimize.c

    r2788 r2908  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-22 05:09:32 $
     11 *  @version $Revision: 1.97 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-01-05 20:59:03 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    190190    n = spline->n;
    191191    klo = p_psVectorBinDisect32(spline->domains, (spline->n)+1, X);
     192    if (klo < 0) {
     193        psLogMsg(__func__, PS_LOG_WARN,
     194                 "WARNING: psMinimize.c: p_psNRSpline1DEval(): p_psVectorBinDisect32 returned an error (%d).\n", klo);
     195        return(NAN);
     196    }
    192197    khi = klo + 1;
    193198    H = (spline->domains)[khi] - (spline->domains)[klo];
     
    245250                                const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
    246251{
    247     PS_PTR_CHECK_NULL(mySpline, NULL);
    248     PS_INT_CHECK_NON_NEGATIVE(mySpline->n, NULL);
    249     PS_VECTOR_CHECK_NULL(x, NULL);
    250252    PS_VECTOR_CHECK_NULL(y, NULL);
    251253    PS_VECTOR_CHECK_TYPE_F32_OR_F64(y, NULL);
     
    278280    // If x==NULL, create an x32 vector with x values set to (0:n).
    279281    if (x == NULL) {
     282
    280283        PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n);
    281284        x32 = x32Static;
     
    286289    PS_VECTOR_CHECK_SIZE_EQUAL(x32, y32, NULL);
    287290    PS_VECTOR_CHECK_SIZE_EQUAL(yErr32, y32, NULL);
     291
     292    /*
     293        XXX:
     294        This can not be implemented until SDR states what order spline should be
     295        created.
     296     
     297        Should we error if mySPline is not NULL?
     298    */
     299    if (mySpline == NULL) {
     300        mySpline = psSpline1DAllocGeneric(x32, 3);
     301    }
     302    PS_PTR_CHECK_NULL(mySpline, NULL);
     303    PS_INT_CHECK_NON_NEGATIVE(mySpline->n, NULL);
     304
    288305
    289306    if (y32->n != (1 + mySpline->n)) {
     
    335352        H = x32->data.F32[i+1] - x32->data.F32[i];
    336353        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    337                 "x data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
     354                "x data (%f - %f) (%f)\n",
     355                x32->data.F32[i],
     356                x32->data.F32[i+1], H);
    338357        //
    339358        // ******** Calculate 0-order term ********
     
    429448        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
    430449    } else {
    431         // XXX: Check size of derivative
     450        PS_IMAGE_CHECK_SIZE(deriv, params->n, coords->n, NULL);
     451        PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    432452    }
    433453
     
    499519    } else {
    500520        PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL);
     521        PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    501522    }
    502523
     
    715736        // Solve A * alpha = Beta
    716737        //
     738        // XXX: How do we know if these functions were successful?
     739        //
    717740        aOut = psMatrixLUD(aOut, perm, A);
    718741        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
     
    963986    }
    964987
     988    // XXX: How do we know if these routines were successful?
    965989    ALUD = psMatrixLUD(ALUD, outPerm, A);
    966990    coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
     
    970994        // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]);
    971995    }
    972 
    973     // for (i=0;i<x->n;i++) {
    974     // printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i],
    975     // myPoly));
    976     // }
    977996
    978997    psFree(A);
     
    10081027
    10091028    psS32 i;
    1010     psPolynomial1D *tmpPoly;
    10111029    psVector *x64 = NULL;
    10121030    psVector *y64 = NULL;
     
    10441062
    10451063    // Call the appropriate vector fitting routine.
     1064    psPolynomial1D *rc = NULL;
    10461065    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    1047         VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     1066        rc = VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    10481067    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    1049         tmpPoly = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
     1068        rc = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    10501069    } else {
    10511070        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    10521071                "unknown polynomial type.\n");
     1072        return(NULL);
     1073    }
     1074    if (rc == NULL) {
     1075        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
    10531076        return(NULL);
    10541077    }
     
    14481471    if (bracket == NULL) {
    14491472        psError(PS_ERR_UNKNOWN, false,
    1450                 "Could not bracket minimum.");
     1473                "Could not bracket minimum.  Returning NAN.\n");
    14511474        return(NAN);
    14521475    }
     
    16361659                if (isnan(mul)) {
    16371660                    psError(PS_ERR_UNKNOWN, false,
    1638                             "Could not perform line minimization");
     1661                            "Could not perform line minimization.  Returning FALSE.\n");
    16391662                    psFree(v);
    16401663                    return(false);
     
    16771700        if (isnan(mul)) {
    16781701            psError(PS_ERR_UNKNOWN, false,
    1679                     "Could not perform line minimization.");
     1702                    "Could not perform line minimization.  Returning FALSE.\n");
    16801703            psFree(v);
    16811704            return(false);
Note: See TracChangeset for help on using the changeset viewer.