IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2235


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

Working with the bracket part of the minimization code.

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r2230 r2235  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-28 23:03:11 $
     11 *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-29 02:00:25 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    11881188}
    11891189
     1190
     1191#define PS_DETERMINE_BRACKET_MAX_ITERATIONS 100
    11901192psVector *p_psDetermineBracket2(psVector *params,
    11911193                                psVector *line,
    … …  
    12001202    float fb = 0.0;
    12011203    float fc = 0.0;
    1202     psS32 iter = 100;
    1203     float aDir = 0.0;
    1204     float cDir = 0.0;
    1205     float new_aDir = 0.0;
    1206     float new_cDir = 0.0;
    1207     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1208     float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
     1204    psS32 iter = 0;
    12091205    psVector *tmp = NULL;
    12101206    psBool boolLineIsNull = true;
    … …  
    12201216        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    12211217                "---- p_psDetermineBracket() end (NULL) ----\n");
    1222         psFree(bracket);
    12231218        return(NULL);
    12241219    }
    12251220
     1221    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    12261222    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
    12271223
    1228     a = 0.0;
    1229     b = stepSize;
    1230 
    1231     PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
    1232     fa = func(tmp, coords);
     1224
     1225    determineINitialBracket(params, paramMask, line, coords, func, &a, &b, &fa, &fb);
     1226
     1227    a = -1.0;
     1228    b = 0.0;
     1229    c = 1.0;
    12331230
    12341231    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
    12351232    fb = func(tmp, coords);
    12361233
    1237     if (fa < fb) {
    1238         aDir = -1;
    1239     } else {
    1240         aDir = 1;
    1241     }
    1242 
     1234    do {
     1235        a*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1236        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
     1237        fa = func(tmp, coords);
     1238    } while (fabs(fb - fa) < FLT_EPSILON);
     1239
     1240    do {
     1241        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1242        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
     1243        fc = func(tmp, coords);
     1244    } while (fabs(fc - fb) < FLT_EPSILON);
     1245
     1246    if ((fb < fa) && (fb < fc)) {
     1247        bracket->data.F32[0] = a;
     1248        bracket->data.F32[1] = b;
     1249        bracket->data.F32[2] = c;
     1250        psFree(tmp);
     1251        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1252                "---- p_psDetermineBracket() end ----\n");
     1253        return(bracket);
     1254    }
     1255
     1256    // The object of this code is to set the initial point a to 0.0, and
     1257    // ensure that (fb < fa).
    12431258    if (fc < fb) {
    1244         cDir = -1;
    1245     } else {
    1246         cDir = 1;
    1247     }
    1248 
    1249     psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1250             "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc);
    1251 
    1252     while (iter > 0) {
     1259        a = b;
     1260        b = c;
     1261        fa = fb;
     1262        fb = fc;
     1263    } else if (fa < fb) {
     1264        float aSave = a;
     1265        float faSave = fa;
     1266
     1267        a = b;
     1268        fa = fb;
     1269
     1270        b = aSave;
     1271        fb = faSave;
     1272    }
     1273    c = b;
     1274
     1275    // At this point we have (a, b) and we know that (fa > fb).  Initially, c=b;
     1276    // We keep stretching b out further from "a" until (fc > previous fc).  If
     1277    // that happens, then we have our bracket.
     1278    iter = 0;
     1279    while (iter < PS_DETERMINE_BRACKET_MAX_ITERATIONS) {
    12531280        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1254                 "psDetermineBracket(): iteration %d\n", iter);
     1281                "psDetermineBracket(): iterationA %d\n", iter);
     1282        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1283
     1284        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
     1285        fc = func(tmp, coords);
     1286
     1287        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1288                "IterA(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1289
    12551290        if ((fb < fa) && (fb < fc)) {
    12561291            bracket->data.F32[0] = a;
    … …  
    12611296                    "---- p_psDetermineBracket() end ----\n");
    12621297            return(bracket);
    1263         }
    1264         stepSize*= (1.0 + stepSize);
    1265         a =- stepSize;
    1266         c =+ stepSize;
    1267 
    1268         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
    1269         fa = func(tmp, coords);
    1270 
    1271         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
    1272         fc = func(tmp, coords);
    1273 
    1274         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1275                 "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    1276 
    1277         if (fa < fb) {
    1278             new_aDir = -1;
    12791298        } else {
    1280             new_aDir = 1;
    1281         }
    1282 
    1283         if (fc < fb) {
    1284             new_cDir = -1;
    1285         } else {
    1286             new_cDir = 1;
    1287         }
    1288         if ((new_aDir == 1) && (aDir == -1)) {
    1289             bracket->data.F32[0] = a;
    1290             bracket->data.F32[1] = b;
    1291             bracket->data.F32[2] = c;
    1292             psFree(tmp);
    1293             psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    1294                     "---- p_psDetermineBracket() end ----\n");
    1295             return(bracket);
    1296         }
    1297 
    1298         if ((new_cDir == 1) && (cDir == -1)) {
    1299             bracket->data.F32[0] = a;
    1300             bracket->data.F32[1] = b;
    1301             bracket->data.F32[2] = c;
    1302             psFree(tmp);
    1303             psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    1304                     "---- p_psDetermineBracket() end ----\n");
    1305             return(bracket);
    1306         }
    1307         aDir = new_aDir;
    1308         cDir = new_cDir;
    1309         iter--;
    1310     }
     1299            b = c;
     1300            fb = fc;
     1301        }
     1302        iter++;
     1303    }
     1304
     1305    // Okay, let's try the other direction.
     1306
    13111307    psFree(tmp);
    13121308    psFree(bracket);
    … …  
    13971393    }
    13981394
    1399     bracket = p_psDetermineBracket(params, line, paramMask, coords, func);
     1395    bracket = p_psDetermineBracket2(params, line, paramMask, coords, func);
    14001396    if (bracket == NULL) {
    14011397        psFree(tmpa);

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

    r2230 r2235  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-28 23:03:11 $
     11 *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-29 02:00:25 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    11881188}
    11891189
     1190
     1191#define PS_DETERMINE_BRACKET_MAX_ITERATIONS 100
    11901192psVector *p_psDetermineBracket2(psVector *params,
    11911193                                psVector *line,
    … …  
    12001202    float fb = 0.0;
    12011203    float fc = 0.0;
    1202     psS32 iter = 100;
    1203     float aDir = 0.0;
    1204     float cDir = 0.0;
    1205     float new_aDir = 0.0;
    1206     float new_cDir = 0.0;
    1207     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1208     float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
     1204    psS32 iter = 0;
    12091205    psVector *tmp = NULL;
    12101206    psBool boolLineIsNull = true;
    … …  
    12201216        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    12211217                "---- p_psDetermineBracket() end (NULL) ----\n");
    1222         psFree(bracket);
    12231218        return(NULL);
    12241219    }
    12251220
     1221    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    12261222    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
    12271223
    1228     a = 0.0;
    1229     b = stepSize;
    1230 
    1231     PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
    1232     fa = func(tmp, coords);
     1224
     1225    determineINitialBracket(params, paramMask, line, coords, func, &a, &b, &fa, &fb);
     1226
     1227    a = -1.0;
     1228    b = 0.0;
     1229    c = 1.0;
    12331230
    12341231    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
    12351232    fb = func(tmp, coords);
    12361233
    1237     if (fa < fb) {
    1238         aDir = -1;
    1239     } else {
    1240         aDir = 1;
    1241     }
    1242 
     1234    do {
     1235        a*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1236        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
     1237        fa = func(tmp, coords);
     1238    } while (fabs(fb - fa) < FLT_EPSILON);
     1239
     1240    do {
     1241        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1242        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
     1243        fc = func(tmp, coords);
     1244    } while (fabs(fc - fb) < FLT_EPSILON);
     1245
     1246    if ((fb < fa) && (fb < fc)) {
     1247        bracket->data.F32[0] = a;
     1248        bracket->data.F32[1] = b;
     1249        bracket->data.F32[2] = c;
     1250        psFree(tmp);
     1251        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1252                "---- p_psDetermineBracket() end ----\n");
     1253        return(bracket);
     1254    }
     1255
     1256    // The object of this code is to set the initial point a to 0.0, and
     1257    // ensure that (fb < fa).
    12431258    if (fc < fb) {
    1244         cDir = -1;
    1245     } else {
    1246         cDir = 1;
    1247     }
    1248 
    1249     psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1250             "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc);
    1251 
    1252     while (iter > 0) {
     1259        a = b;
     1260        b = c;
     1261        fa = fb;
     1262        fb = fc;
     1263    } else if (fa < fb) {
     1264        float aSave = a;
     1265        float faSave = fa;
     1266
     1267        a = b;
     1268        fa = fb;
     1269
     1270        b = aSave;
     1271        fb = faSave;
     1272    }
     1273    c = b;
     1274
     1275    // At this point we have (a, b) and we know that (fa > fb).  Initially, c=b;
     1276    // We keep stretching b out further from "a" until (fc > previous fc).  If
     1277    // that happens, then we have our bracket.
     1278    iter = 0;
     1279    while (iter < PS_DETERMINE_BRACKET_MAX_ITERATIONS) {
    12531280        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1254                 "psDetermineBracket(): iteration %d\n", iter);
     1281                "psDetermineBracket(): iterationA %d\n", iter);
     1282        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
     1283
     1284        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
     1285        fc = func(tmp, coords);
     1286
     1287        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1288                "IterA(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1289
    12551290        if ((fb < fa) && (fb < fc)) {
    12561291            bracket->data.F32[0] = a;
    … …  
    12611296                    "---- p_psDetermineBracket() end ----\n");
    12621297            return(bracket);
    1263         }
    1264         stepSize*= (1.0 + stepSize);
    1265         a =- stepSize;
    1266         c =+ stepSize;
    1267 
    1268         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
    1269         fa = func(tmp, coords);
    1270 
    1271         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
    1272         fc = func(tmp, coords);
    1273 
    1274         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    1275                 "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    1276 
    1277         if (fa < fb) {
    1278             new_aDir = -1;
    12791298        } else {
    1280             new_aDir = 1;
    1281         }
    1282 
    1283         if (fc < fb) {
    1284             new_cDir = -1;
    1285         } else {
    1286             new_cDir = 1;
    1287         }
    1288         if ((new_aDir == 1) && (aDir == -1)) {
    1289             bracket->data.F32[0] = a;
    1290             bracket->data.F32[1] = b;
    1291             bracket->data.F32[2] = c;
    1292             psFree(tmp);
    1293             psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    1294                     "---- p_psDetermineBracket() end ----\n");
    1295             return(bracket);
    1296         }
    1297 
    1298         if ((new_cDir == 1) && (cDir == -1)) {
    1299             bracket->data.F32[0] = a;
    1300             bracket->data.F32[1] = b;
    1301             bracket->data.F32[2] = c;
    1302             psFree(tmp);
    1303             psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
    1304                     "---- p_psDetermineBracket() end ----\n");
    1305             return(bracket);
    1306         }
    1307         aDir = new_aDir;
    1308         cDir = new_cDir;
    1309         iter--;
    1310     }
     1299            b = c;
     1300            fb = fc;
     1301        }
     1302        iter++;
     1303    }
     1304
     1305    // Okay, let's try the other direction.
     1306
    13111307    psFree(tmp);
    13121308    psFree(bracket);
    … …  
    13971393    }
    13981394
    1399     bracket = p_psDetermineBracket(params, line, paramMask, coords, func);
     1395    bracket = p_psDetermineBracket2(params, line, paramMask, coords, func);
    14001396    if (bracket == NULL) {
    14011397        psFree(tmpa);
Note: See TracChangeset for help on using the changeset viewer.