IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

Modifying the minimization algorithms.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2221 r2228  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.72 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-27 23:31:43 $
     11 *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-28 22:46:57 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    408408}
    409409
     410
     411/******************************************************************************
     412XXX: We assume unnormalized gaussians.
     413XXX: Currently, yErr is ignored.
     414 *****************************************************************************/
     415psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
     416                                  const psVector *params,
     417                                  const psArray *coords)
     418{
     419    psTrace(".psLib.dataManip.psMinimize", 4,
     420            "---- psMinimizeLMChi2Gauss1D() begin ----\n");
     421
     422    PS_PTR_CHECK_NULL(coords, NULL);
     423    PS_PTR_CHECK_NULL(params, NULL);
     424    float x;
     425    int i;
     426    float mean = params->data.F32[0];
     427    float stdev = params->data.F32[1];
     428    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
     429
     430    psTrace(".psLib.dataManip.psMinimize", 6,
     431            "(mean, stdev) is (%f, %f)\n", mean, stdev);
     432
     433    if (deriv == NULL) {
     434        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
     435    }
     436
     437    for (i=0;i<coords->n;i++) {
     438        x = ((psVector *) (coords->data[i]))->data.F32[0];
     439        out->data.F32[i] = psGaussian(x, mean, stdev, false);
     440        //        psTrace(".psLib.dataManip.psMinimize", 6,
     441        //            "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]);
     442    }
     443
     444    for (i=0;i<coords->n;i++) {
     445        x = ((psVector *) (coords->data[i]))->data.F32[0];
     446        float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
     447        deriv->data.F32[i][0] = tmp / (stdev * stdev);
     448
     449        tmp = (x - mean) * (x - mean) *
     450              psGaussian(x, mean, stdev, 0);
     451        deriv->data.F32[i][1] = tmp / (stdev * stdev * stdev);
     452    }
     453
     454    psTrace(".psLib.dataManip.psMinimize", 4,
     455            "---- psMinimizeLMChi2Gauss1D() end ----\n");
     456    return(out);
     457}
     458
     459
    410460/******************************************************************************
    411461psMinimizeLMChi2():  This routine will take an procedure which calculates
    … …  
    429479XXX: This must work for both F32 and F64.  F32 is currently implemented.
    430480     Note: since the LUD routines are only implemented in F64, then we
    431      will have to convert all F32 input vectors to F64 regardless.  So,
     481>     will have to convert all F32 input vectors to F64 regardless.  So,
    432482     the F64 port might be.
    433483 *****************************************************************************/
    … …  
    457507    }
    458508
    459     psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     509    psTrace(".psLib.dataManip.psMinimize", 4,
    460510            "---- psMinimizeLMChi2() begin ----\n");
    461511    psS32 numData = y->n;
    … …  
    490540    float newChi2 = 0.0;
    491541    float lamda = 0.00005;
    492 
    493     psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     542    lamda = 0.05;
     543
     544    psTrace(".psLib.dataManip.psMinimize", 6,
    494545            "min->maxIter is %d\n", min->maxIter);
    495     psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     546    psTrace(".psLib.dataManip.psMinimize", 6,
    496547            "min->tol is %f\n", min->tol);
    497548
    … …  
    504555
    505556    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
    506         psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     557        psTrace(".psLib.dataManip.psMinimize", 4,
    507558                "------------------------------------------------------\n");
    508         psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     559        psTrace(".psLib.dataManip.psMinimize", 4,
    509560                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
    510561
    … …  
    515566        currValueVec = func(deriv, params, x);
    516567        for (n=0;n<numData;n++) {
    517             currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
     568            currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) *
     569                       (currValueVec->data.F32[n] - y->data.F32[n]);
     570            psTrace(".psLib.dataManip.psMinimize", 6,
     571                    "data[%d], chi2 calculation+= (%f * %f)\n", n,
     572                    currValueVec->data.F32[n], y->data.F32[n]);
    518573        }
    519574
    520575        for (p=0;p<numParams;p++) {
    521             psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     576            psTrace(".psLib.dataManip.psMinimize", 6,
    522577                    "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
    523578        }
    524         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     579        psTrace(".psLib.dataManip.psMinimize", 6,
    525580                "Current chi-squared is (%f)\n", currChi2);
    526581
    … …  
    548603            // XXX: multiple by -1 here?
    549604            (beta->data.F64[p])*= -1.0;
    550             psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     605            psTrace(".psLib.dataManip.psMinimize", 6,
    551606                    "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
    552607        }
    … …  
    580635        }
    581636        for (j=0;j<numParams;j++) {
    582             psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
     637            psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n");
    583638            for (k=0;k<numParams;k++) {
    584                 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);
    585             }
    586             psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
     639                psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]);
     640            }
     641            psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n");
    587642        }
    588643
    … …  
    597652        //
    598653        for (i=0;i<numParams;i++) {
    599             psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     654            psTrace(".psLib.dataManip.psMinimize", 6,
    600655                    "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
    601656            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    … …  
    607662        }
    608663
    609         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     664        psTrace(".psLib.dataManip.psMinimize", 6,
    610665                "Calling func() with new parameters:\n");
    611666        for (i=0;i<numParams;i++) {
    612             psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     667            psTrace(".psLib.dataManip.psMinimize", 6,
    613668                    "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
    614669        }
    … …  
    621676        newValueVec = func(deriv, newParams, x);
    622677        for (n=0;n<numData;n++) {
    623             newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
     678            newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) *
     679                      (newValueVec->data.F32[n] - y->data.F32[n]);
     680
    624681        }
    625682        psFree(newValueVec);
    626683
    627         psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     684        psTrace(".psLib.dataManip.psMinimize", 4,
    628685                "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
    629686
    … …  
    640697            }
    641698            lamda*= 0.1;
     699            psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n");
    642700        } else {
    643701            lamda*= 10.0;
    644         }
    645         psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     702            psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n");
     703        }
     704        psTrace(".psLib.dataManip.psMinimize", 4,
    646705                "lamda is %f\n", lamda);
    647706        min->iter++;
    … …  
    662721    }
    663722
    664     psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     723    psTrace(".psLib.dataManip.psMinimize", 4,
    665724            "---- psMinimizeLMChi2() end (false) ----\n");
    666725    return(false);
    … …  
    9421001}
    9431002
     1003#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
     1004for (int i=0;i<BASE->n;i++) { \
     1005    if (BASEMASK->data.U8[i] == 0) { \
     1006        OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \
     1007    } else { \
     1008        OUT->data.F32[i] = BASE->data.F32[i]; \
     1009    } \
     1010} \
     1011
     1012
    9441013/******************************************************************************
    9451014p_psDetermineBracket():  This routine takes as input an arbitrary function,
    … …  
    9471016function produces as output a bracket [a, b, c] such that
    9481017 
    949         f(param + b * line) is less than f(param + a * line) and
    950                                          f(param + c * line).
    951  
    952 Algorithm: XXX completely ad hoc: start with the user-supplied starting
     1018f(param + b * line) is less than f(param + a * line) and
     1019f(param + c * line).
     1020 
     1021Algorithm:
     1022XXX completely ad hoc:
     1023start with the user-supplied starting
    9531024parameter and call that b.  Calculate a/c as a fractional amount
    9541025smaller/larger than b.  Repeat this process until a local minimum is found.
    9551026 
    956 XXX: new algorithm: start at x=0, expand in one direction until the function
     1027XXX:
     1028new algorithm:
     1029start at x=0, expand in one direction until the function
    9571030decreases.  Then you have two points in the bracket.  Keep going until it
    9581031increases, or x is too large.  If thst does not work, expand in the other
    9591032direction.
    9601033 
    961 XXX: This is F32 only
    962  *****************************************************************************/
     1034XXX:
     1035This is F32 only.
     1036 
     1037XXX:
     1038output bracket vector should be an input as well.
     1039*****************************************************************************/
    9631040psVector *p_psDetermineBracket(psVector *params,
    9641041                               psVector *line,
    … …  
    9661043                               const psArray *coords,
    9671044                               psMinimizePowellFunc func)
     1045{
     1046    float a = 0.0;
     1047    float b = 0.0;
     1048    float c = 0.0;
     1049    float fa = 0.0;
     1050    float fb = 0.0;
     1051    float fc = 0.0;
     1052    psS32 iter = 100;
     1053    float aDir = 0.0;
     1054    float cDir = 0.0;
     1055    float new_aDir = 0.0;
     1056    float new_cDir = 0.0;
     1057    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
     1058    float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
     1059    psVector *tmp = NULL;
     1060    psS32 i = 0;
     1061    psS32 null = 0;
     1062
     1063    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1064            "---- p_psDetermineBracket() begin ----\n");
     1065
     1066    for (i=0;i<params->n;i++) {
     1067        if (paramMask->data.U8[i] == 0) {
     1068            if (line->data.F32[i] >= FLT_EPSILON) {
     1069                null = 1;
     1070            }
     1071        }
     1072    }
     1073
     1074    if (null == 0) {
     1075        psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
     1076                "p_psDetermineBracket() called with zero line vector.\n");
     1077        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1078                "---- p_psDetermineBracket() end (NULL) ----\n");
     1079        psFree(bracket);
     1080        return(NULL);
     1081    }
     1082
     1083    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
     1084
     1085    b = 0;
     1086    a = -stepSize;
     1087    c = stepSize;
     1088
     1089    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
     1090    fa = func(tmp, coords);
     1091
     1092    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
     1093    fb = func(tmp, coords);
     1094
     1095    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
     1096    fc = func(tmp, coords);
     1097
     1098    if (fa < fb) {
     1099        aDir = -1;
     1100    } else {
     1101        aDir = 1;
     1102    }
     1103
     1104    if (fc < fb) {
     1105        cDir = -1;
     1106    } else {
     1107        cDir = 1;
     1108    }
     1109
     1110    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1111            "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc);
     1112
     1113    while (iter > 0) {
     1114        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1115                "psDetermineBracket(): iteration %d\n", iter);
     1116        if ((fb < fa) && (fb < fc)) {
     1117            bracket->data.F32[0] = a;
     1118            bracket->data.F32[1] = b;
     1119            bracket->data.F32[2] = c;
     1120            psFree(tmp);
     1121            psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1122                    "---- p_psDetermineBracket() end ----\n");
     1123            return(bracket);
     1124        }
     1125        stepSize*= (1.0 + stepSize);
     1126        a =- stepSize;
     1127        c =+ stepSize;
     1128
     1129        for (i=0;i<params->n;i++) {
     1130            if (paramMask->data.U8[i] == 0) {
     1131                tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1132            } else {
     1133                tmp->data.F32[i] = params->data.F32[i];
     1134            }
     1135        }
     1136        fa = func(tmp, coords);
     1137
     1138        for (i=0;i<params->n;i++) {
     1139            if (paramMask->data.U8[i] == 0) {
     1140                tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1141            } else {
     1142                tmp->data.F32[i] = params->data.F32[i];
     1143            }
     1144        }
     1145        fc = func(tmp, coords);
     1146
     1147        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1148                "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1149
     1150        if (fa < fb) {
     1151            new_aDir = -1;
     1152        } else {
     1153            new_aDir = 1;
     1154        }
     1155
     1156        if (fc < fb) {
     1157            new_cDir = -1;
     1158        } else {
     1159            new_cDir = 1;
     1160        }
     1161        if ((new_aDir == 1) && (aDir == -1)) {
     1162            bracket->data.F32[0] = a;
     1163            bracket->data.F32[1] = b;
     1164            bracket->data.F32[2] = c;
     1165            psFree(tmp);
     1166            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1167                    "---- p_psDetermineBracket() end ----\n");
     1168            return(bracket);
     1169        }
     1170
     1171        if ((new_cDir == 1) && (cDir == -1)) {
     1172            bracket->data.F32[0] = a;
     1173            bracket->data.F32[1] = b;
     1174            bracket->data.F32[2] = c;
     1175            psFree(tmp);
     1176            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1177                    "---- p_psDetermineBracket() end ----\n");
     1178            return(bracket);
     1179        }
     1180        aDir = new_aDir;
     1181        cDir = new_cDir;
     1182        iter--;
     1183    }
     1184    psFree(tmp);
     1185    psFree(bracket);
     1186    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1187            "---- p_psDetermineBracket() end (NULL) ----\n");
     1188    return(NULL);
     1189}
     1190
     1191psVector *p_psDetermineBracket2(psVector *params,
     1192                                psVector *line,
     1193                                const psVector *paramMask,
     1194                                const psArray *coords,
     1195                                psMinimizePowellFunc func)
    9681196{
    9691197    float a = 0.0;
    … …  
    13651593    float fqp = 0.0;
    13661594    float diff = 0.0;
     1595    int iterationNumber = 0;
    13671596
    13681597    psTrace(".psLib.dataManip.psMinimizePowell", 4,
    … …  
    14071636    }
    14081637
    1409     min->iter = 0;
    1410     while (min->iter < min->maxIter) {
    1411         min->iter++;
     1638    while (iterationNumber < min->maxIter) {
     1639        iterationNumber++;
    14121640        psTrace(".psLib.dataManip.psMinimizePowell", 6,
    1413                 "psMinimizePowell() iteration %d\n", min->iter);
     1641                "psMinimizePowell() iteration %d\n", iterationNumber);
    14141642
    14151643        // 3: For each dimension in params, move Q only in the vector v[i] to
    … …  
    14181646        baseFuncVal = func(Q, coords);
    14191647        currFuncVal = baseFuncVal;
     1648        psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1649                "Current function value is %f\n", currFuncVal);
     1650
    14201651        diff = 0.0;
    14211652        biggestDiff = 0;
    … …  
    14231654        for (i=0;i<numDims;i++) {
    14241655            if (myParamMask->data.U8[i] == 0) {
    1425                 dummyMin.maxIter = 100;
     1656                dummyMin.maxIter = 20;
    14261657                dummyMin.tol = 0.01;
    14271658                mul = p_psLineMin(&dummyMin, Q, v[i], myParamMask, coords, func);
    … …  
    14331664            }
    14341665        }
    1435 
     1666        psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1667                "New function value is %f\n", currFuncVal);
    14361668        // 4: Set the vector u = Q - P
    14371669        for (i=0;i<numDims;i++) {
    … …  
    14501682
    14511683        // 5: Move Q only in the direction u, and minimize the function.
    1452         psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
    14531684        for (i=0;i<numDims;i++) {
    14541685            psTrace(".psLib.dataManip.psMinimizePowell", 6,
    … …  
    14701701                psFree(myParamMask);
    14711702            }
     1703            min->iter = iterationNumber;
    14721704            psTrace(".psLib.dataManip.psMinimizePowell", 4,
    14731705                    "---- psMinimizePowell() end (true) ----\n");
    … …  
    15131745            psFree(pQP);
    15141746            psFree(Q);
     1747            if (paramMask == NULL) {
     1748                psFree(myParamMask);
     1749            }
     1750            min->iter = iterationNumber;
    15151751            psTrace(".psLib.dataManip.psMinimizePowell", 4,
    15161752                    "---- psMinimizePowell() end (true) ----\n");
    1517             if (paramMask == NULL) {
    1518                 psFree(myParamMask);
    1519             }
    15201753            return(true);
    15211754        }
    … …  
    15321765    psFree(pQP);
    15331766    psFree(Q);
     1767    min->iter = iterationNumber;
    15341768    psTrace(".psLib.dataManip.psMinimizePowell", 4,
    15351769            "---- psMinimizePowell() end (false) ----\n");
    … …  
    15491783                       const psArray *coords)
    15501784{
     1785    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
     1786            "---- myPowellChi2Func() begin ----\n");
     1787    PS_VECTOR_CHECK_NULL(params, NAN);
     1788    PS_VECTOR_CHECK_EMPTY(params, NAN);
     1789    PS_VECTOR_CHECK_NULL(myValue, NAN);
     1790    PS_VECTOR_CHECK_EMPTY(myValue, NAN);
     1791    PS_PTR_CHECK_NULL(coords, NAN);
     1792
    15511793    float chi2 = 0.0;
    15521794    float d;
    … …  
    15551797
    15561798    tmp = Chi2PowellFunc(params, coords);
    1557     for (i=0;i<coords->n;i++) {
    1558         d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i];
    1559         chi2+= d * d;
     1799    if (myError == NULL) {
     1800        for (i=0;i<coords->n;i++) {
     1801            d = (tmp->data.F32[i] - myValue->data.F32[i]);
     1802            chi2+= d * d;
     1803        }
     1804    } else {
     1805        for (i=0;i<coords->n;i++) {
     1806            d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i];
     1807            chi2+= d * d;
     1808        }
    15601809    }
    15611810    psFree(tmp);
     1811    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
     1812            "---- myPowellChi2Func() end (%f) ----\n", chi2);
    15621813    return(chi2);
    15631814}
    … …  
    15841835    myValue = (psVector *) value;
    15851836    myError = (psVector *) error;
     1837
    15861838    Chi2PowellFunc = func;
    15871839
Note: See TracChangeset for help on using the changeset viewer.