IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 21, 2004, 7:09:32 PM (22 years ago)
Author:
gusciora
Message:

Changed types to PS types.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2741 r2788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.95 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-17 00:18:31 $
     11 *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-22 05:09:32 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 *  XXX: must follow coding name standards on local functions.
    17  *
    18  *  XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via
    19  *       LMM) is not addressed here.  We are waiting for subsequent SDRs
    20  *       which will redefine the LMM functions.
    2117 *
    2218 */
    … …  
    6460XXX: Use a static vector.
    6561 *****************************************************************************/
    66 void psBuildSums1D(double x,
     62void psBuildSums1D(psF64 x,
    6763                   psS32 polyOrder,
    6864                   psVector* sums)
    6965{
    7066    psS32 i = 0;
    71     double xSum = 0.0;
     67    psF64 xSum = 0.0;
    7268
    7369    if (sums == NULL) {
    … …  
    10096XXX: do an F64 version?
    10197 *****************************************************************************/
    102 float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
     98psF32 *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    10399                             const psVector* restrict y)        ///< Coordinates
    104100{
    … …  
    108104    psS32 i;
    109105    psS32 k;
    110     float sig;
    111     float p;
     106    psF32 sig;
     107    psF32 p;
    112108    psS32 n = y->n;
    113     float *u = (float *) psAlloc(n * sizeof(float));
    114     float *derivs2 = (float *) psAlloc(n * sizeof(float));
    115     float *X = (float *) & (x->data.F32[0]);
    116     float *Y = (float *) & (y->data.F32[0]);
    117     float qn;
     109    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
     110    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
     111    psF32 *X = (psF32 *) & (x->data.F32[0]);
     112    psF32 *Y = (psF32 *) & (y->data.F32[0]);
     113    psF32 qn;
    118114
    119115    // XXX: The second derivatives at the endpoints, undefined in the SDR,
    … …  
    168164 *****************************************************************************/
    169165/*
    170 float p_psNRSpline1DEval(psSpline1D *spline,
     166psF32 p_psNRSpline1DEval(psSpline1D *spline,
    171167                         const psVector* restrict x,
    172168                         const psVector* restrict y,
    173                          float X)
     169                         psF32 X)
    174170{
    175171    PS_PTR_CHECK_NULL(spline, NAN);
    … …  
    185181    psS32 klo;
    186182    psS32 khi;
    187     float H;
    188     float A;
    189     float B;
    190     float C;
    191     float D;
    192     float Y;
     183    psF32 H;
     184    psF32 A;
     185    psF32 B;
     186    psF32 C;
     187    psF32 D;
     188    psF32 Y;
    193189 
    194190    n = spline->n;
    … …  
    258254            "---- psVectorFitSpline1D() begin ----\n");
    259255    psS32 numSplines = (y->n)-1;
    260     float tmp;
    261     float H;
     256    psF32 tmp;
     257    psF32 H;
    262258    psS32 i;
    263     float slope;
     259    psF32 slope;
    264260    psVector *x32 = NULL;
    265261    psVector *y32 = NULL;
    … …  
    421417    psTrace(".psLib.dataManip.psMinimize", 4,
    422418            "---- psMinimizeLMChi2Gauss1D() begin ----\n");
    423     float x;
    424     int i;
    425     float mean = params->data.F32[0];
    426     float stdev = params->data.F32[1];
     419    psF32 x;
     420    psS32 i;
     421    psF32 mean = params->data.F32[0];
     422    psF32 stdev = params->data.F32[1];
    427423    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    428424
    … …  
    443439    for (i=0;i<coords->n;i++) {
    444440        x = ((psVector *) (coords->data[i]))->data.F32[0];
    445         float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
     441        psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false);
    446442        deriv->data.F32[i][0] = tmp / (stdev * stdev);
    447443        tmp = (x - mean) * (x - mean) *
    … …  
    491487    PS_PTR_CHECK_NULL(params, NULL);
    492488
    493     double normalization = params->data.F32[0];
    494     double x0 = params->data.F32[1];
    495     double y0 = params->data.F32[2];
    496     double sigmaX = params->data.F32[3];
    497     double sigmaY = params->data.F32[4];
    498     double theta = params->data.F32[5];
     489    psF64 normalization = params->data.F32[0];
     490    psF64 x0 = params->data.F32[1];
     491    psF64 y0 = params->data.F32[2];
     492    psF64 sigmaX = params->data.F32[3];
     493    psF64 sigmaY = params->data.F32[4];
     494    psF64 theta = params->data.F32[5];
    499495    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    500496
    … …  
    502498        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
    503499    } else {
    504         // XXX: Check size of derivative
     500        PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL);
    505501    }
    506502
    … …  
    508504            "---- psMinimizeLMChi2Gauss2D() begin ----\n");
    509505
    510     for (int i=0;i<coords->n;i++) {
    511         double x = ((psVector *) coords->data[i])->data.F32[0];
    512         double y = ((psVector *) coords->data[i])->data.F32[0];
    513 
    514         double u = - (x-x0)*cos(theta) + (y-y0)*sin(theta);
    515         double v = (x-x0)*cos(theta) + (y-y0)*sin(theta);
    516 
    517         double flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +
    518                                              v*v/(2.0 * sigmaY * sigmaY)))/
    519                       (2.0 * M_PI * sigmaX * sigmaY);
     506    for (psS32 i=0;i<coords->n;i++) {
     507        psF64 x = ((psVector *) coords->data[i])->data.F32[0];
     508        psF64 y = ((psVector *) coords->data[i])->data.F32[0];
     509
     510        psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta);
     511        psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta);
     512
     513        psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +
     514                                            v*v/(2.0 * sigmaY * sigmaY)))/
     515                     (2.0 * M_PI * sigmaX * sigmaY);
    520516        out->data.F32[i] = flux;
    521517
    … …  
    584580        PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL);
    585581    }
    586 
    587     // XXX: Generate code that modifies covar matrix if not-NULL.
    588582    if (covar != NULL) {
    589583        PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL);
    … …  
    610604    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    611605    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    612 
    613     //    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
    614     //    for (i=0;i<numData;i++) {
    615     //        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
    616     //    }
    617606    psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32);
    618 
    619     psVector *currValueVec;
    620     psVector *newValueVec;
    621 
    622     float currChi2 = 0.0;
    623     float newChi2 = 0.0;
    624     float lamda = 0.00005;
     607    psVector *currValueVec = NULL;
     608    psVector *newValueVec = NULL;
     609    psF32 currChi2 = 0.0;
     610    psF32 newChi2 = 0.0;
     611    psF32 lamda = 0.00005;
    625612    lamda = 0.05;
    626613
    … …  
    711698                if (j == k) {
    712699                    A->data.F64[j][k] =
    713                         (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
     700                        (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]);
    714701                } else {
    715                     A->data.F64[j][k] = (double) alpha->data.F32[j][k];
     702                    A->data.F64[j][k] = (psF64) alpha->data.F32[j][k];
    716703                }
    717704            }
    … …  
    741728            } else {
    742729                newParams->data.F32[i] = params->data.F32[i] -
    743                                          (float) paramDeltasF64->data.F64[i];
     730                                         (psF32) paramDeltasF64->data.F64[i];
    744731            }
    745732        }
    … …  
    777764            // We already masked params.
    778765            for (i=0;i<numParams;i++) {
    779                 params->data.F32[i] = (float) newParams->data.F32[i];
     766                params->data.F32[i] = (psF32) newParams->data.F32[i];
    780767            }
    781768            lamda*= 0.1;
    … …  
    826813    psS32 k;
    827814    psS32 n = x->n;
    828     double fac;
    829     double sum;
     815    psF64 fac;
     816    psF64 sum;
    830817    PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64);
    831818    psScalar *fScalar;
    … …  
    836823    // variable declarations.  I retain them here to maintain coherence
    837824    // with the NR code.
    838     double min = -1.0;
    839     double max = 1.0;
    840     double bma = 0.5 * (max-min);  // 1
    841     double bpa = 0.5 * (max+min);  // 0
     825    psF64 min = -1.0;
     826    psF64 max = 1.0;
     827    psF64 bma = 0.5 * (max-min);  // 1
     828    psF64 bpa = 0.5 * (max+min);  // 0
    842829
    843830    // In this loop, we first calculate the values of X for which the
    … …  
    851838    for (psS32 i=0;i<n;i++) {
    852839        // NR 5.8.4
    853         double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
    854         double X = (Y + bma + bpa) - 1.0;
     840        psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n));
     841        psF64 X = (Y + bma + bpa) - 1.0;
    855842        tmpScalar.data.F64 = X;
    856843
    … …  
    873860    // coefficients of the Chebyshev polynomial: NR 5.8.7.
    874861
    875     fac = 2.0/((float) n);
     862    fac = 2.0/((psF32) n);
    876863    // XXX: is this loop bound correct?
    877864    for (j=0;j<myPoly->n;j++) {
    … …  
    879866        for (k=0;k<n;k++) {
    880867            sum+= f->data.F64[k] *
    881                   cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n));
     868                  cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n));
    882869        }
    883870
    … …  
    10751062 *****************************************************************************/
    10761063psMinimization *psMinimizationAlloc(psS32 maxIter,
    1077                                     float tol)
     1064                                    psF32 tol)
    10781065{
    10791066    PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL);
    … …  
    10921079// LINE to it.  We assume BASEMASK is non-null.
    10931080#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
    1094 for (int i=0;i<BASE->n;i++) { \
     1081for (psS32 i=0;i<BASE->n;i++) { \
    10951082    if (BASEMASK->data.U8[i] == 0) { \
    10961083        OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \
    … …  
    11021089#define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \
    11031090BOOL_VAR = true; \
    1104 for (int i=0;i<IN->n;i++) { \
     1091for (psS32 i=0;i<IN->n;i++) { \
    11051092    if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \
    11061093        BOOL_VAR = false; \
    … …  
    11111098#define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \
    11121099BOOL_VAR = true; \
    1113 for (int i=0;i<IN->n;i++) { \
     1100for (psS32 i=0;i<IN->n;i++) { \
    11141101    if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \
    11151102        BOOL_VAR = false; \
    … …  
    11521139                               psMinimizePowellFunc func)
    11531140{
    1154     float a = 0.0;
    1155     float b = 0.0;
    1156     float c = 0.0;
    1157     float fa = 0.0;
    1158     float fb = 0.0;
    1159     float fc = 0.0;
     1141    psF32 a = 0.0;
     1142    psF32 b = 0.0;
     1143    psF32 c = 0.0;
     1144    psF32 fa = 0.0;
     1145    psF32 fb = 0.0;
     1146    psF32 fc = 0.0;
    11601147    psS32 iter = 100;
    1161     float aDir = 0.0;
    1162     float cDir = 0.0;
    1163     float new_aDir = 0.0;
    1164     float new_cDir = 0.0;
     1148    psF32 aDir = 0.0;
     1149    psF32 cDir = 0.0;
     1150    psF32 new_aDir = 0.0;
     1151    psF32 new_cDir = 0.0;
    11651152    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1166     float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
     1153    psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
    11671154    psVector *tmp = NULL;
    11681155    psBool boolLineIsNull = true;
    … …  
    13011288                                psMinimizePowellFunc func)
    13021289{
    1303     float a = 0.0;
    1304     float b = 0.0;
    1305     float c = 0.0;
    1306     float fa = 0.0;
    1307     float fb = 0.0;
    1308     float fc = 0.0;
     1290    psF32 a = 0.0;
     1291    psF32 b = 0.0;
     1292    psF32 c = 0.0;
     1293    psF32 fa = 0.0;
     1294    psF32 fb = 0.0;
     1295    psF32 fc = 0.0;
    13091296    psS32 iter = 0;
    13101297    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
    13111298    psBool boolLineIsNull = true;
    1312     float prevMin = 0.0;
    1313     int countMin = 0;
     1299    psF32 prevMin = 0.0;
     1300    psS32 countMin = 0;
    13141301
    13151302    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    … …  
    14031390 *****************************************************************************/
    14041391#define PS_LINEMIN_MAX_ITERATIONS 30
    1405 float p_psLineMin(psMinimization *min,
     1392psF32 p_psLineMin(psMinimization *min,
    14061393                  psVector *params,
    14071394                  psVector *line,
    … …  
    14231410    PS_PTR_CHECK_NULL(func, NAN);
    14241411    psVector *bracket;
    1425     float a = 0.0;
    1426     float b = 0.0;
    1427     float c = 0.0;
    1428     float n = 0.0;
    1429     float fa = 0.0;
    1430     float fb = 0.0;
    1431     float fc = 0.0;
    1432     float fn = 0.0;
    1433     float mul = 0.0;
     1412    psF32 a = 0.0;
     1413    psF32 b = 0.0;
     1414    psF32 c = 0.0;
     1415    psF32 n = 0.0;
     1416    psF32 fa = 0.0;
     1417    psF32 fb = 0.0;
     1418    psF32 fc = 0.0;
     1419    psF32 fn = 0.0;
     1420    psF32 mul = 0.0;
    14341421    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
    14351422    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
    … …  
    15441531This routine must minimize a possibly multi-dimensional function.  The
    15451532function to be minimized "func" is:
    1546     float func(psVector *params, psArray *coords)
     1533    psF32 func(psVector *params, psArray *coords)
    15471534The "params" are the parameters of the function which are varied.  The data
    15481535points at which the function is varied are in the argument "coords" which is
    … …  
    15781565    psVector *myParamMask = NULL;
    15791566    psMinimization dummyMin;
    1580     float mul = 0.0;
    1581     float baseFuncVal = 0.0;
    1582     float currFuncVal = 0.0;
     1567    psF32 mul = 0.0;
     1568    psF32 baseFuncVal = 0.0;
     1569    psF32 currFuncVal = 0.0;
    15831570    psS32 biggestIter = 0;
    1584     float biggestDiff = 0.0;
    1585     int iterationNumber = 0;
     1571    psF32 biggestDiff = 0.0;
     1572    psS32 iterationNumber = 0;
    15861573
    15871574    psTrace(".psLib.dataManip.psMinimizePowell", 4,
    … …  
    17111698            }
    17121699        }
    1713         float fqp = func(pQP, coords);
    1714         float term1 = (baseFuncVal - currFuncVal) - biggestDiff;
     1700        psF32 fqp = func(pQP, coords);
     1701        psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff;
    17151702        term1*= term1;
    17161703        term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp);
    1717         float term2 = baseFuncVal - fqp;
     1704        psF32 term2 = baseFuncVal - fqp;
    17181705        term2*= term2 * biggestDiff;
    17191706        if (term1 < term2) {
    … …  
    17601747    PS_PTR_CHECK_NULL(params, NULL);
    17611748
    1762     float x;
    1763     int i;
    1764     float mean = params->data.F32[0];
    1765     float stdev = params->data.F32[1];
     1749    psF32 x;
     1750    psS32 i;
     1751    psF32 mean = params->data.F32[0];
     1752    psF32 stdev = params->data.F32[1];
    17661753    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    17671754
    … …  
    17831770XXX: This is F32 only
    17841771 *****************************************************************************/
    1785 float myPowellChi2Func(const psVector *params,
     1772psF32 myPowellChi2Func(const psVector *params,
    17861773                       const psArray *coords)
    17871774{
    … …  
    17941781    PS_PTR_CHECK_NULL(coords, NAN);
    17951782
    1796     float chi2 = 0.0;
    1797     float d;
     1783    psF32 chi2 = 0.0;
     1784    psF32 d;
    17981785    psS32 i;
    17991786    psVector *tmp;
Note: See TracChangeset for help on using the changeset viewer.