IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 1, 2004, 1:57:08 PM (22 years ago)
Author:
gusciora
Message:

Fixed the robust stats routines for stdev/mean. Fixed minimization routines
for fitting gaussians. CVS:

Enter Log. Lines beginning with `CVS:' are removed automatically CVS: CVS:
Committing in . CVS: CVS: Modified Files: CVS:
src/dataManip/psMinimize.c src/dataManip/psMinimize.h CVS:
src/dataManip/psStats.c src/dataManip/makedir/psMinimize.d CVS:
src/dataManip/makedir/psStats.d CVS: test/dataManip/tst_psMinimize07.c

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2247 r2250  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-11-01 18:57:05 $
     11 *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-11-01 23:57:08 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    572572                       (currValueVec->data.F32[n] - y->data.F32[n]);
    573573            psTrace(".psLib.dataManip.psMinimize", 6,
    574                     "data[%d], chi2 calculation+= (%f * %f)\n", n,
     574                    "data[%d], chi2 calculation+= (%f - %f)^2\n", n,
    575575                    currValueVec->data.F32[n], y->data.F32[n]);
    576576        }
    … …  
    10041004}
    10051005
     1006// This macro takes as input the vector BASE and adds a multiple of the vector
     1007// LINE to it.  We assume BASEMASK is non-null.
    10061008#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
    10071009for (int i=0;i<BASE->n;i++) { \
    … …  
    10361038and the parameter to vary, and the line along which it must vary.  This
    10371039function produces as output a bracket [a, b, c] such that
    1038  
    1039 f(param + b * line) is less than f(param + a * line) and
    1040 f(param + c * line).
     1040f(param + b * line) < f(param + a * line)
     1041f(param + b * line) < f(param + c * line)
     1042a < b < c
    10411043 
    10421044Algorithm:
    … …  
    11921194
    11931195
     1196#define RETURN_FINAL_BRACKET(d) \
     1197if (a < c) { \
     1198    bracket->data.F32[0] = a; \
     1199    bracket->data.F32[1] = b; \
     1200    bracket->data.F32[2] = c; \
     1201} else { \
     1202    bracket->data.F32[0] = c; \
     1203    bracket->data.F32[1] = b; \
     1204    bracket->data.F32[2] = a; \
     1205} \
     1206psTrace(".psLib.dataManip.p_psDetermineBracket", 4, \
     1207        "---- p_psDetermineBracket() end ----\n"); \
     1208psTrace(".psLib.dataManip.p_psDetermineBracket", 4, "Final bracket (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); \
     1209return(bracket); \
     1210
    11941211#define PS_DETERMINE_BRACKET_MAX_ITERATIONS 100
    11951212psVector *p_psDetermineBracket2(psVector *params,
    … …  
    12061223    float fc = 0.0;
    12071224    psS32 iter = 0;
    1208     psVector *tmp = NULL;
     1225    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
    12091226    psBool boolLineIsNull = true;
    12101227    float prevMin = 0.0;
    12111228    int countMin = 0;
    1212 
    1213     psTrace(".PS_SEG.PS_PWD.PS_FILE.__func__", 6,
    1214             "---- __func__() begin COOL");
    1215 
    1216     psTrace(".psLib.dataManip.__func__", 4,
    1217             "---- p_psDetermineBracket() COOL ----\n");
    12181229
    12191230    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    … …  
    12301241    }
    12311242
    1232     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1233     tmp = psVectorAlloc(params->n, PS_TYPE_F32);
    1234     a = -0.5;
    1235     b = 0.0;
    1236     c = 0.5;
    1237 
    1238     PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
    1239     fb = func(tmp, coords);
    1240 
     1243    // We determine in what x-direction does the function decrease.
     1244    a = 0.0;
     1245    fa = func(params, coords);
     1246    b = 0.5;
     1247    iter = 0;
    12411248    do {
    1242         a*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
    1243         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
    1244         fa = func(tmp, coords);
    1245     } while (fabs(fb - fa) < FLT_EPSILON);
    1246 
    1247     do {
    1248         c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
    1249         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
    1250         fc = func(tmp, coords);
    1251     } while (fabs(fc - fb) < FLT_EPSILON);
    1252 
    1253     psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1254             "initial (fa, Fb, Fc) is (%f, %f, %f)\n", fa, fb, fc);
    1255 
    1256     if ((fb < fa) && (fb < fc)) {
    1257         bracket->data.F32[0] = a;
    1258         bracket->data.F32[1] = b;
    1259         bracket->data.F32[2] = c;
    1260         psFree(tmp);
    1261         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1262                 "---- p_psDetermineBracket() end ----\n");
    1263         return(bracket);
    1264     }
    1265 
    1266     // The object of this code is to set the initial point a to 0.0, and
    1267     // ensure that (fb < fa).
    1268     if (fc < fb) {
    1269         a = b;
    1270         b = c;
    1271         fa = fb;
    1272         fb = fc;
    1273     } else if (fa < fb) {
    1274         float aSave = a;
    1275         float faSave = fa;
    1276 
     1249        b*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1250        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
     1251        fb = func(tmp, coords);
     1252    } while ((fabs(fb - fa) < FLT_EPSILON) && (iter++ < 100));
     1253
     1254    if (fb > fa) {
    12771255        a = b;
    12781256        fa = fb;
    1279 
    1280         b = aSave;
    1281         fb = faSave;
     1257        b = 0.0;
     1258        fb = func(params, coords);
    12821259    }
    12831260    c = b;
    12841261
    1285     // At this point we have (a, b) and we know that (fa > fb).  Initially, c=b;
     1262    // At this point we have (a, b) and we know that (fa >= fb).  Initially, c=b;
    12861263    // We keep stretching b out further from "a" until (fc > previous fc).  If
    12871264    // that happens, then we have our bracket.
     1265    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    12881266    iter = 0;
    12891267    while (iter < PS_DETERMINE_BRACKET_MAX_ITERATIONS) {
    12901268        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    12911269                "psDetermineBracket(): iterationA %d\n", iter);
    1292         c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1270        c+= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE) * (c - a);
    12931271
    12941272        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
    … …  
    12961274
    12971275        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1298                 "IterA(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1276                "Iteration(%d) (bracket): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    12991277
    13001278        if ((fb < fa) && (fb < fc)) {
    1301             bracket->data.F32[0] = a;
    1302             bracket->data.F32[1] = b;
    1303             bracket->data.F32[2] = c;
    1304             psFree(tmp);
    1305             psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1306                     "---- p_psDetermineBracket() end ----\n");
    1307             return(bracket);
     1279            RETURN_FINAL_BRACKET();
    13081280        } else {
    13091281            b = c;
    … …  
    13121284
    13131285        // This code maintains a count of how many times the minimum fc has
    1314         // statyed the same.  If it gets too high, we exit this loop.
     1286        // stayed the same.  If it gets too high, we exit this loop.
    13151287        if (fc == prevMin) {
    13161288            countMin++;
    … …  
    13201292        prevMin = fc;
    13211293        if (countMin == 10) {
    1322             bracket->data.F32[0] = a;
    1323             bracket->data.F32[1] = b;
    1324             bracket->data.F32[2] = c;
    1325             psFree(tmp);
    1326             psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1327                     "---- p_psDetermineBracket() end ----\n");
    1328             return(bracket);
     1294            RETURN_FINAL_BRACKET();
    13291295        }
    13301296
    … …  
    13321298    }
    13331299
    1334     // Okay, let's try the other direction.
    1335 
    1336     psFree(tmp);
    13371300    psFree(bracket);
    13381301    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    1339             "---- p_psDetermineBracket() end (NULL) ----\n");
     1302            "---- p_psDetermineBracket() end (NULL) (BAD) ----\n");
    13401303    return(NULL);
    13411304}
    … …  
    13511314function.
    13521315XXX: This is F32 only
     1316XXX: Since this is an internal function, many of the parameter checks are
     1317     redundant.
    13531318 *****************************************************************************/
     1319#define PS_LINEMIN_MAX_ITERATIONS 30
    13541320float p_psLineMin(psMinimization *min,
    13551321                  psVector *params,
    … …  
    13591325                  psMinimizePowellFunc func)
    13601326{
     1327    PS_PTR_CHECK_NULL(min, NAN);
     1328    PS_VECTOR_CHECK_NULL(params, NAN);
     1329    PS_VECTOR_CHECK_EMPTY(params, NAN);
     1330    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1331    PS_VECTOR_CHECK_NULL(line, NAN);
     1332    PS_VECTOR_CHECK_EMPTY(line, NAN);
     1333    PS_VECTOR_CHECK_TYPE(line, PS_TYPE_F32, NAN);
     1334    PS_VECTOR_CHECK_NULL(paramMask, NAN);
     1335    PS_VECTOR_CHECK_EMPTY(paramMask, NAN);
     1336    PS_VECTOR_CHECK_TYPE(paramMask, PS_TYPE_U8, NAN);
     1337    PS_PTR_CHECK_NULL(coords, NAN);
     1338    PS_PTR_CHECK_NULL(func, NAN);
    13611339    psVector *bracket;
    13621340    float a = 0.0;
    … …  
    13691347    float fn = 0.0;
    13701348    float mul = 0.0;
    1371     psVector *tmpa = NULL;
    1372     psVector *tmpb = NULL;
    1373     psVector *tmpc = NULL;
    1374     psVector *tmpn = NULL;
     1349    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
     1350    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
     1351    PS_VECTOR_GEN_STATIC_RECYCLED(tmpc, params->n, PS_TYPE_F32);
     1352    PS_VECTOR_GEN_STATIC_RECYCLED(tmpn, params->n, PS_TYPE_F32);
    13751353    psS32 i = 0;
    13761354    psS32 boolLineIsNull = true;
    1377 
    1378     psTrace(".psLib.dataManip.p_psLineMin", 4,
    1379             "---- p_psLineMin() begin ----\n");
    1380     psTrace(".psLib.dataManip.p_psLineMin", 6,
    1381             "min->maxIter is %d\n", min->maxIter);
    1382     psTrace(".psLib.dataManip.p_psLineMin", 6,
    1383             "min->tol is %f\n", min->tol);
    1384 
    1385     if (paramMask != NULL) {
    1386         for (i=0;i<params->n;i++) {
    1387             if (paramMask->data.U8[i] == 0) {
    1388                 if (line->data.F32[i] >= FLT_EPSILON) {
    1389                     boolLineIsNull = false;
    1390                     psTrace(".psLib.dataManip.p_psLineMin", 4,
    1391                             "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
    1392                 }
    1393             }
    1394         }
    1395     } else {
    1396         for (i=0;i<params->n;i++) {
    1397             if (line->data.F32[i] >= FLT_EPSILON) {
    1398                 boolLineIsNull = false;
    1399                 psTrace(".psLib.dataManip.p_psLineMin", 4,
    1400                         "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
    1401             }
    1402         }
    1403     }
     1355    psS32 numIterations = 0;
     1356
     1357    psTrace(".psLib.dataManip.p_psLineMin", 4, "---- p_psLineMin() begin ----\n");
     1358    PS_VECTOR_F32_CHECK_ZERO_VECTOR(line, boolLineIsNull);
    14041359
    14051360    if (boolLineIsNull == true) {
     1361        min->value = func(params, coords);
    14061362        psTrace(".psLib.dataManip.p_psLineMin", 2,
    1407                 "p_psLineMin() called with zero line vector.\n");
     1363                "p_psLineMin() called with zero line vector.  Return 0.0.  Function valus is %f\n", min->value);
    14081364        return(0.0);
    14091365    }
    1410 
    1411     tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
    1412     tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
    1413     tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
    1414     tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
    14151366
    14161367    for (i=0;i<params->n;i++) {
    14171368        psTrace(".psLib.dataManip.p_psLineMin", 6,
    1418                 "params, paramMask, line [%d] is (%f %d %f)\n", i,
     1369                "(params, paramMask, line)[%d] is (%f %d %f)\n", i,
    14191370                params->data.F32[i],
    14201371                paramMask->data.U8[i],
    … …  
    14241375    bracket = p_psDetermineBracket2(params, line, paramMask, coords, func);
    14251376    if (bracket == NULL) {
    1426         psFree(tmpa);
    1427         psFree(tmpb);
    1428         psFree(tmpc);
    1429         psFree(tmpn);
    14301377        psError(__func__, "(1) Could not bracket minimum.");
    14311378        return(NAN);
    14321379    }
    1433 
    1434     min->iter = 0;
    1435     while (min->iter < min->maxIter) {
    1436         min->iter++;
     1380    numIterations = 0;
     1381    while (numIterations < PS_LINEMIN_MAX_ITERATIONS) {
     1382        numIterations++;
    14371383        psTrace(".psLib.dataManip.p_psLineMin", 6,
    1438                 "p_psLineMin(): iteration %d\n", min->iter);
     1384                "p_psLineMin(): iteration %d\n", numIterations);
    14391385
    14401386        a = bracket->data.F32[0];
    14411387        b = bracket->data.F32[1];
    14421388        c = bracket->data.F32[2];
    1443 
    1444         for (i=0;i<params->n;i++) {
    1445             if (paramMask->data.U8[i] == 0) {
    1446                 tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
    1447                 tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
    1448                 tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
    1449             } else {
    1450                 tmpa->data.F32[i] = params->data.F32[i];
    1451                 tmpb->data.F32[i] = params->data.F32[i];
    1452                 tmpc->data.F32[i] = params->data.F32[i];
    1453             }
    1454         }
     1389        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpa, a);
     1390        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpb, b);
     1391        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpc, c);
    14551392        fa = func(tmpa, coords);
    14561393        fb = func(tmpb, coords);
    14571394        fc = func(tmpc, coords);
    14581395        psTrace(".psLib.dataManip.p_psLineMin", 6,
    1459                 "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter,
    1460                 a, b, c, fa, fb, fc);
     1396                "LineMin: f(%f %f %f) is (%f %f %f)\n", a, b, c, fa, fb, fc);
    14611397
    14621398        // We determine which is the biggest segment in [a,b,c] then split
    … …  
    14951431        }
    14961432        psTrace(".psLib.dataManip.p_psLineMin", 6,
    1497                 "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
     1433                "LineMin: new bracket is (%f %f %f)\n", bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
    14981434
    14991435        mul = bracket->data.F32[1];
    15001436        if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
    1501             for (i=0;i<params->n;i++) {
    1502                 if (paramMask->data.U8[i] == 0) {
    1503                     params->data.F32[i]+= mul * line->data.F32[i];
    1504                 }
    1505             }
     1437            PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
    15061438            min->value = func(params, coords);
    15071439            psFree(bracket);
    1508             psFree(tmpa);
    1509             psFree(tmpb);
    1510             psFree(tmpc);
    1511             psFree(tmpn);
    15121440            psTrace(".psLib.dataManip.p_psLineMin", 4,
    1513                     "---- p_psLineMin() end.a (%f) ----\n", mul);
     1441                    "---- p_psLineMin() end.a (%f) (%f) ----\n", mul, min->value);
    15141442            return(mul);
    15151443        }
    15161444    }
    15171445
     1446    mul = bracket->data.F32[1];
     1447    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
     1448    min->value = func(params, coords);
     1449    psTrace(".psLib.dataManip.p_psLineMin", 4,
     1450            "---- p_psLineMin() end.b (%f) %f ----\n", mul, min->value);
     1451
    15181452    psFree(bracket);
    1519     psFree(tmpa);
    1520     psFree(tmpb);
    1521     psFree(tmpc);
    1522     psFree(tmpn);
    1523 
    1524     psTrace(".psLib.dataManip.p_psLineMin", 4,
    1525             "---- p_psLineMin() end.b (0.0) ----\n");
    1526     return(0.0);
     1453    return(mul);
    15271454}
    15281455
    … …  
    15601487
    15611488    psS32 numDims = params->n;
     1489    PS_VECTOR_GEN_STATIC_RECYCLED(pQP, numDims, PS_TYPE_F32);
     1490    PS_VECTOR_GEN_STATIC_RECYCLED(u, numDims, PS_TYPE_F32);
     1491    PS_VECTOR_GEN_STATIC_RECYCLED(Q, numDims, PS_TYPE_F32);
    15621492    psS32 i = 0;
    15631493    psS32 j = 0;
    15641494    psVector **v = NULL;
    1565     psVector *pQP = NULL;
    1566     psVector *Q = NULL;
    1567     psVector *u = NULL;
    15681495    psVector *myParamMask = NULL;  // XXX: Use a recylced vector here.
    15691496    psMinimization dummyMin;
    … …  
    15861513            "min->tol is %f\n", min->tol);
    15871514
    1588     pQP = psVectorAlloc(numDims, PS_TYPE_F32);
    1589     u = psVectorAlloc(numDims, PS_TYPE_F32);
    1590     Q = psVectorAlloc(numDims, PS_TYPE_F32);
    1591 
    15921515    if (paramMask == NULL) {
    1593         myParamMask = psVectorAlloc(params->n, PS_TYPE_U8);
     1516        myParamMask = psVectorRecycle(myParamMask, params->n, PS_TYPE_U8);
     1517        p_psMemSetPersistent(myParamMask, true);
     1518        p_psMemSetPersistent(myParamMask->data.V, true);
    15941519        for (i=0;i<myParamMask->n;i++) {
    15951520            myParamMask->data.U8[i] = 0;
    … …  
    16431568                if (isnan(mul)) {
    16441569                    psError(__func__, "Could not perform line minimization (1).\n");
    1645                     psFree(pQP);
    1646                     psFree(u);
    1647                     psFree(Q);
    1648                     if (paramMask == NULL) {
    1649                         psFree(myParamMask);
    1650                     }
    1651                     for (i=0;i<numDims;i++) {
    1652                         psFree(v[i]);
     1570                    for (j=0;j<numDims;j++) {
     1571                        psFree(v[j]);
    16531572                    }
    16541573                    psFree(v);
    16551574                    return(false);
    16561575                }
     1576                psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1577                        "LineMin along dimension %d has multiple %f\n", i, mul);
     1578
    16571579                if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
    16581580                    biggestDiff = fabs(dummyMin.value - currFuncVal);
    … …  
    16641586        psTrace(".psLib.dataManip.psMinimizePowell", 6,
    16651587                "New function value is %f\n", currFuncVal);
     1588
    16661589        // 4: Set the vector u = Q - P
    16671590        for (i=0;i<numDims;i++) {
    … …  
    16701593
    16711594                psTrace(".psLib.dataManip.psMinimizePowell", 6,
    1672                         "u=Q-P (%f = %f - %f)\n", u->data.F32[i],
     1595                        "u[i]=Q[i]-P[i] (%f = %f - %f)\n", u->data.F32[i],
    16731596                        Q->data.F32[i],
    16741597                        params->data.F32[i]);
    … …  
    16881611        if (isnan(mul)) {
    16891612            psError(__func__, "Could not perform line minimization. (2)\n");
    1690             psFree(pQP);
    1691             psFree(u);
    1692             psFree(Q);
    1693             if (paramMask == NULL) {
    1694                 psFree(myParamMask);
    1695             }
    16961613            for (i=0;i<numDims;i++) {
    16971614                psFree(v[i]);
    … …  
    17071624            }
    17081625            psFree(v);
    1709             psFree(u);
    1710             psFree(pQP);
    1711             psFree(Q);
    1712             if (paramMask == NULL) {
    1713                 psFree(myParamMask);
    1714             }
    17151626            min->iter = iterationNumber;
    17161627            psTrace(".psLib.dataManip.psMinimizePowell", 4,
    1717                     "---- psMinimizePowell() end (true) ----\n");
     1628                    "---- psMinimizePowell() end (1)(true) ----\n");
    17181629            return(true);
    17191630        }
    … …  
    17541665            }
    17551666            psFree(v);
    1756             psFree(u);
    1757             psFree(pQP);
    1758             psFree(Q);
    1759             if (paramMask == NULL) {
    1760                 psFree(myParamMask);
    1761             }
    17621667            min->iter = iterationNumber;
    17631668            psTrace(".psLib.dataManip.psMinimizePowell", 4,
    1764                     "---- psMinimizePowell() end (true) ----\n");
     1669                    "---- psMinimizePowell() end (2) (true) ----\n");
    17651670            return(true);
    17661671        }
    17671672    }
    17681673
    1769     if (paramMask == NULL) {
    1770         psFree(myParamMask);
    1771     }
    17721674    for (i=0;i<numDims;i++) {
    17731675        psFree(v[i]);
    17741676    }
    17751677    psFree(v);
    1776     psFree(u);
    1777     psFree(pQP);
    1778     psFree(Q);
    17791678    min->iter = iterationNumber;
    17801679    psTrace(".psLib.dataManip.psMinimizePowell", 4,
    1781             "---- psMinimizePowell() end (false) ----\n");
     1680            "---- psMinimizePowell() end (0) (false) ----\n");
    17821681    return(false);
     1682}
     1683
     1684
     1685/******************************************************************************
     1686XXX: We assume unnormalized gaussians.
     1687XXX: Currently, yErr is ignored.
     1688 *****************************************************************************/
     1689psVector *psMinimizePowellChi2Gauss1D(const psVector *params,
     1690                                      const psArray *coords)
     1691{
     1692    PS_PTR_CHECK_NULL(coords, NULL);
     1693    PS_PTR_CHECK_NULL(params, NULL);
     1694
     1695    float x;
     1696    int i;
     1697    float mean = params->data.F32[0];
     1698    float stdev = params->data.F32[1];
     1699    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
     1700
     1701    for (i=0;i<coords->n;i++) {
     1702        x = ((psVector *) (coords->data[i]))->data.F32[0];
     1703        out->data.F32[i] = psGaussian(x, mean, stdev, false);
     1704    }
     1705
     1706    return(out);
    17831707}
    17841708
    … …  
    18221746    psFree(tmp);
    18231747    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
    1824             "---- myPowellChi2Func() end (%f) ----\n", chi2);
     1748            "---- myPowellChi2Func() end (chi2 is %f) ----\n", chi2);
    18251749    return(chi2);
    18261750}
Note: See TracChangeset for help on using the changeset viewer.