IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5813


Ignore:
Timestamp:
Dec 19, 2005, 1:58:47 PM (20 years ago)
Author:
gusciora
Message:

Implemented the chi-squared algorithm for calculated Chebyshev poly
coefficients.

Location:
trunk/psLib
Files:
3 added
4 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psCoord.c

    r5682 r5813  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-12-05 21:33:29 $
     12*  @version $Revision: 1.97 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-12-19 23:58:47 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6161should rename this.
    6262 
    63 XXX: Use the ADD version which is based on determinants.
    64  
    6563XXX: This code no longer makes sense.  The merge must be reviewed.
    6664 *****************************************************************************/
    … …  
    8583
    8684    if (0) {
    87         //This is sample code from IfA.  It didn't work initially, and I did not
    88         //spend any time debugging it.
    89         psF64 a = transform->x->coeff[1][0];
    90         psF64 b = transform->x->coeff[0][1];
    91         psF64 c = transform->y->coeff[1][0];
    92         psF64 d = transform->y->coeff[0][1];
    93         psF64 e = transform->x->coeff[0][0];
    94         psF64 f = transform->y->coeff[0][0];
    95 
    96         psF64 invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
    97 
    98         // Not entirely sure why this works, but it appears to do so
    99         out->x->coeff[1][0] = invDet * a;
    100         out->x->coeff[0][1] = - invDet * b;
    101         out->y->coeff[1][0] = - invDet * c;
    102         out->y->coeff[0][1] = invDet * d;
    103 
    104         out->x->coeff[0][0] = - invDet * (d * e + c * f);
    105         out->y->coeff[0][0] = - invDet * (b * e + a * f);
    106     }
    107 
    108     if (1) {
    109         // XXX: There is no reason to execute this code and the following code.
     85        // XXX: Get rid of this code.
    11086        psF64 A = transform->x->coeff[1][0];
    11187        psF64 B = transform->x->coeff[0][1];
    … …  
    129105    }
    130106
    131 
    132     // XXX: There is no reason to execute this code and the previous codes.
    133     // unless the cross terms are available, set these matrix elements to 0
    134107    psF64 r12 = 0.0;
    135108    if (transform->x->nY == 1) {

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

    r5624 r5813  
    77*  polynomials.  It also contains a Gaussian functions.
    88*
    9 *  @version $Revision: 1.133 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-11-30 02:00:09 $
     9*  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-12-19 23:58:47 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    165165not nOrder.
    166166 *****************************************************************************/
    167 static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly)
     167psPolynomial1D **createChebyshevPolys(psS32 numPolys)
     168{
     169    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(numPolys, 1, NULL);
     170
     171    if (0) {
     172        psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
     173        for (psS32 i = 0; i < numPolys; i++) {
     174            chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
     175            chebPolys[i]->coeff[i] = 1;
     176        }
     177        return (chebPolys);
     178    } else {
     179        psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
     180        for (psS32 i = 0; i < numPolys; i++) {
     181            chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
     182        }
     183
     184        // Create the Chebyshev polynomials.
     185        // Polynomial i has i-th order.
     186        chebPolys[0]->coeff[0] = 1.0;
     187        if (numPolys >= 2) {
     188            chebPolys[1]->coeff[1] = 1.0;
     189
     190            for (psS32 i = 2; i < numPolys; i++) {
     191                for (psS32 j = 0; j < chebPolys[i - 1]->nX+1; j++) {
     192                    chebPolys[i]->coeff[j + 1] = 2.0 * chebPolys[i - 1]->coeff[j];
     193                }
     194                for (psS32 j = 0; j < chebPolys[i - 2]->nX+1; j++) {
     195                    chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
     196                }
     197            }
     198        }
     199
     200        return (chebPolys);
     201    }
     202}
     203
     204/*
     205static psPolynomial1D **createChebyshevPolysOld(psS32 maxChebyPoly)
    168206{
    169207    PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL);
    170 
     208 
    171209    psPolynomial1D **chebPolys = NULL;
    172 
     210 
    173211    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
    174212    for (psS32 i = 0; i < maxChebyPoly; i++) {
    175213        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
    176214    }
    177 
     215 
    178216    // Create the Chebyshev polynomials.
    179217    // Polynomial i has i-th order.
    180218    chebPolys[0]->coeff[0] = 1;
    181 
     219 
    182220    // XXX: Bug 296
    183221    if (maxChebyPoly > 1) {
    184222        chebPolys[1]->coeff[1] = 1;
    185 
     223 
    186224        for (psS32 i = 2; i < maxChebyPoly; i++) {
    187225            for (psS32 j = 0; j < chebPolys[i - 1]->nX; j++) {
    … …  
    196234        printf("WARNING: %u-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
    197235    }
    198 
     236 
    199237    return (chebPolys);
    200238}
     239*/
    201240
    202241/*****************************************************************************
    … …  
    236275// XXX: How does the mask vector effect Crenshaw's formula?
    237276// XXX: We assume that x is scaled between -1.0 and 1.0;
    238 static psF64 chebPolynomial1DEval(psF64 x,
    239                                   const psPolynomial1D* poly)
    240 {
    241     PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    242     // XXX: Create a macro for this in psConstants.h
    243     if (poly->nX < 1) {
    244         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial is order %u.", poly->nX);
    245         return(NAN);
    246     }
     277// XXX: Create a faster version for low-order Chebyshevs.
     278static psF64 chebPolynomial1DEval(
     279    psF64 x,
     280    const psPolynomial1D* poly)
     281{
     282    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, NAN);
     283    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(poly->nX, 0, NAN);
    247284    psVector *d;
     285
    248286    // XXX: n should be nTerms here (for clarity).  Or get rid of the variable.
    249     unsigned int n = 1 + poly->nX;
     287    unsigned int nTerms = 1 + poly->nX;
    250288    unsigned int i;
    251289    psF64 tmp = 0.0;
    252290
    253291    // Special case where the Chebyshev poly is constant.
    254     if (n == 1) {
     292    if (nTerms == 1) {
    255293        if (poly->mask[0] == 0) {
    256294            tmp += poly->coeff[0];
    … …  
    260298
    261299    // Special case where the Chebyshev poly is linear.
    262     if (n == 2) {
     300    if (nTerms == 2) {
    263301        if (poly->mask[0] == 0) {
    264302            tmp+= poly->coeff[0];
    … …  
    270308    }
    271309
    272     // General case where the Chebyshev poly has 2 or more terms.
    273     d = psVectorAlloc(n, PS_TYPE_F64);
    274     if(poly->mask[n-1] == 0) {
    275         d->data.F64[n-1] = poly->coeff[n-1];
     310    if (1) {
     311        // General case where the Chebyshev poly has 2 or more terms.
     312        d = psVectorAlloc(nTerms, PS_TYPE_F64);
     313        if(poly->mask[nTerms-1] == 0) {
     314            d->data.F64[nTerms-1] = poly->coeff[nTerms-1];
     315        } else {
     316            d->data.F64[nTerms-1] = 0.0;
     317        }
     318
     319        d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]);
     320        if(poly->mask[nTerms-2] == 0) {
     321            d->data.F64[nTerms-2] += poly->coeff[nTerms-2];
     322        }
     323
     324        for (i=nTerms-3;i>=1;i--) {
     325            d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
     326                             (d->data.F64[i+2]);
     327            if(poly->mask[i] == 0) {
     328                d->data.F64[i] += poly->coeff[i];
     329            }
     330        }
     331
     332        tmp = (x * d->data.F64[1]) -
     333              (d->data.F64[2]);
     334        if(poly->mask[0] == 0) {
     335            tmp += (0.5 * poly->coeff[0]);
     336        }
     337        psFree(d);
    276338    } else {
    277         d->data.F64[n-1] = 0.0;
    278     }
    279 
    280     d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]);
    281     if(poly->mask[n-2] == 0) {
    282         d->data.F64[n-2] += poly->coeff[n-2];
    283     }
    284 
    285     for (i=n-3;i>=1;i--) {
    286         d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
    287                          (d->data.F64[i+2]);
    288         if(poly->mask[i] == 0) {
    289             d->data.F64[i] += poly->coeff[i];
    290         }
    291     }
    292 
    293     tmp = (x * d->data.F64[1]) -
    294           (d->data.F64[2]);
    295     if(poly->mask[0] == 0) {
    296         tmp += (0.5 * poly->coeff[0]);
    297     }
    298     psFree(d);
     339        // This is old code that does not use Clenshaw's formula.  Get rid of it.
     340        psPolynomial1D **chebPolys = createChebyshevPolys(1 + poly->nX);
     341
     342        tmp = 0.0;
     343        for (psS32 i=0;i<(1 + poly->nX);i++) {
     344            tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x));
     345        }
     346        tmp-= (poly->coeff[0]/2.0);
     347
     348        for (psS32 i=0;i<(1 + poly->nX);i++) {
     349            psFree(chebPolys[i]);
     350        }
     351        psFree(chebPolys);
     352    }
     353
    299354    return(tmp);
    300 
    301     /* This is old code that does not use Clenshaw's formula.  Get rid of it.
    302 
    303     psS32 i;
    304     psF32 tmp;
    305     psPolynomial1D **chebPolys = NULL;
    306 
    307     chebPolys = createChebyshevPolys(1 + poly->nX);
    308 
    309     tmp = 0.0;
    310     for (i=0;i<(1 + poly->nX);i++) {
    311         tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i]));
    312     }
    313     tmp-= (poly->coeff[0]/2.0);
    314 
    315 
    316     return(tmp);
    317     */
    318355}
    319356
    … …  
    628665    newPoly->type = type;
    629666    newPoly->nX = n;
     667    newPoly->p_min = NAN;
     668    newPoly->p_max = NAN;
    630669    newPoly->coeff = psAlloc((n + 1) * sizeof(psF64));
    631670    newPoly->coeffErr = psAlloc((n + 1) * sizeof(psF64));

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

    r5294 r5813  
    1111 *  @author GLG, MHPCC
    1212 *
    13  *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-10-12 21:02:20 $
     13 *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-12-19 23:58:47 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    7777    psPolynomialType type;              ///< Polynomial type
    7878    unsigned int nX;             ///< Polynomial order
     79    psF64 p_min;
     80    psF64 p_max;
    7981    psF64 *coeff;                       ///< Coefficients
    8082    psF64 *coeffErr;                    ///< Error in coefficients
    … …  
    298300);
    299301
     302
     303// XXX: Coding Standard
     304psPolynomial1D **createChebyshevPolys(psS32 numPolys);
     305
     306typedef struct
     307{
     308    int n;                              ///< The number of Chebyshev polys.
     309    psPolynomial1D **chebyPolys;        ///< THe chebyshev polys
     310
     311}
     312p_chebyPolys;
     313
     314
     315
    300316/** \} */ // End of MathGroup Functions
    301317

  • trunk/psLib/test/math/Makefile.am

    r5155 r5813  
    2626    tst_psMatrixVectorArithmetic03 \
    2727    tst_psMatrixVectorArithmetic04 \
    28     tst_psMinimize04 \
    29     tst_psMinimize04_F32 \
    30     tst_psMinimize04b \
    31     tst_psMinimize04b_F32 \
    3228    tst_psMinimize05 \
    3329    tst_psMinimize06 \
    … …  
    4440    tst_psStats09 \
    4541    tst_psSpline1D \
    46     tst_psRandom
     42    tst_psRandom \
     43    tst_psPolyFit1DCheby \
     44    tst_psPolyFit1DOrd \
     45    tst_psPolyFit2DOrd
     46
    4747
    4848tst_psFunc00_SOURCES =  tst_psFunc00.c
    … …  
    6767tst_psMatrixVectorArithmetic03_SOURCES =  tst_psMatrixVectorArithmetic03.c
    6868tst_psMatrixVectorArithmetic04_SOURCES =  tst_psMatrixVectorArithmetic04.c
    69 tst_psMinimize04_SOURCES =  tst_psMinimize04.c
    70 tst_psMinimize04_F32_SOURCES =  tst_psMinimize04_F32.c
    71 tst_psMinimize04b_SOURCES =  tst_psMinimize04b.c
    72 tst_psMinimize04b_F32_SOURCES =  tst_psMinimize04b_F32.c
    7369tst_psMinimizeVector2D_F64_SOURCES = tst_psMinimizeVector2D_F64.c
    7470tst_psMinimizeVector2D_F32_SOURCES = tst_psMinimizeVector2D_F32.c
    … …  
    8581tst_psStats09_SOURCES =  tst_psStats09.c
    8682tst_psRandom_SOURCES =  tst_psRandom.c
     83tst_psPolyFit1DCheby_SOURCES = tst_psPolyFit1DCheby.c
     84tst_psPolyFit1DOrd_SOURCES = tst_psPolyFit1DOrd.c
     85tst_psPolyFit2DOrd_SOURCES = tst_psPolyFit2DOrd.c
    8786
    8887check_DATA =

  • trunk/psLib/test/math/tst_psStats07.c

    r5155 r5813  
    2020psS32 t00()
    2121{
     22
    2223    psStats * myStats = NULL;
    2324    psS32 testStatus = true;
    … …  
    510511{
    511512    psLogSetFormat("HLNM");
     513    psLogSetLevel(PS_LOG_INFO);
    512514    //
    513515    // We list pertinent psStats.c functions here for debugging ease.
Note: See TracChangeset for help on using the changeset viewer.