IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1719


Ignore:
Timestamp:
Sep 7, 2004, 8:03:26 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r1718 r1719  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 06:00:49 $
     9 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-08 06:01:21 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1313 */
    14 
    1514/*****************************************************************************/
    16 
    1715/*  INCLUDE FILES                                                            */
    18 
    1916/*****************************************************************************/
    2017#include <stdlib.h>
    … …  
    3734
    3835/*****************************************************************************/
    39 
    4036/* DEFINE STATEMENTS                                                         */
    41 
    4237/*****************************************************************************/
    4338
    … …  
    8378
    8479/*****************************************************************************/
    85 
    8680/* TYPE DEFINITIONS                                                          */
    87 
    8881/*****************************************************************************/
    8982static void polynomial1DFree(psPolynomial1D* myPoly);
    … …  
    9790
    9891/*****************************************************************************/
    99 
    10092/* GLOBAL VARIABLES                                                          */
    101 
    10293/*****************************************************************************/
    10394
    … …  
    10596
    10697/*****************************************************************************/
    107 
    10898/* FILE STATIC VARIABLES                                                     */
    109 
    11099/*****************************************************************************/
    111100
    … …  
    113102
    114103/*****************************************************************************/
    115 
    116104/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    117 
    118105/*****************************************************************************/
    119106
    … …  
    121108
    122109/*****************************************************************************/
    123 
    124110/*  FUNCTION IMPLEMENTATION - PUBLIC                                         */
    125 
    126111/*****************************************************************************/
    127112

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

    r1629 r1719  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-08-26 20:10:25 $
     11 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-08 06:02:23 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4343#include "psMinimize.h"
    4444#include "psMatrix.h"
     45#include "psTrace.h"
    4546
    4647/*****************************************************************************/
    … …  
    569570    gsl_vector *x;
    570571
     572    psTrace(".psLib.dataManip.psMinimize", 4,
     573            "---- psMinimize() begin ----\n");
     574
    571575    PS_CHECK_NULL_VECTOR(initialGuess);
    572576    PS_CHECK_EMPTY_VECTOR(initialGuess);
    … …  
    655659        }
    656660    }
     661    psTrace(".psLib.dataManip.psMinimize", 4,
     662            "---- psMinimize() end ----\n");
    657663    return (initialGuess);
    658664}
    … …  
    691697    float chiSqOld = 0.0;
    692698
     699    psTrace(".psLib.dataManip.psMinimizeChi2", 4,
     700            "---- psMinimizeChi2() begin ----\n");
     701
    693702    PS_CHECK_NULL_IMAGE(domain);
    694703    PS_CHECK_EMPTY_IMAGE(domain);
    … …  
    837846
    838847    // Bye bye.
     848    psTrace(".psLib.dataManip.psMinimizeChi2", 4,
     849            "---- psMinimizeChi2() end ----\n");
    839850    return (initialGuess);
    840851}
    … …  
    863874    psVector* xSums = NULL;
    864875
     876    psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
     877            "---- psVectorFitPolynomial1D() begin ----\n");
    865878    // printf("psVectorFitPolynomial1D()\n");
    866879    // for (i=0;i<x->n;i++) {
    … …  
    936949    psFree(xSums);
    937950
     951    psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
     952            "---- psVectorFitPolynomial1D() begin ----\n");
    938953    return (myPoly);
     954}
     955
     956psMinimization *psMinimizationAlloc(int maxIter,
     957                                    float tol)
     958{
     959    psMinimization *min = psAlloc(sizeof(psMinimization));
     960    min->maxIter = maxIter;
     961    min->tol = tol;
     962    min->value = 0.0;
     963    min->iter = 0;
     964    min->lastDelta = 0.0;
     965
     966    return(min);
    939967}
    940968
    941969#define STEP_SIZE 0.10
    942970/******************************************************************************
    943     This routine takes as input an arbitrary function of a single parameter.
    944     This function produces as output a bracket [a, b, c] such that f(b) is
    945     less than f(a) and f(b).
     971    This routine takes as input an arbitrary function, and the parameter to
     972    vary.  This function produces as output a bracket [a, b, c] such that
     973    f(b) is less than f(a) and f(b).
    946974 
    947975    Algorithm: XXX completely ad hoc: start with the user-supplied starting
    … …  
    950978    found.
    951979 *****************************************************************************/
    952 psVector *p_psDetermineBracket(psVector *params,
    953                                const psArray *coords,
    954                                psMinimizePowellFunc func)
     980psVector *p_psDetermineBracketOld(psVector *params,
     981                                  int dim,
     982                                  const psArray *coords,
     983                                  psMinimizePowellFunc func)
    955984{
    956985    float a = 0.0;
    … …  
    966995    float new_cDir = 0.0;
    967996    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    968     psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
    969 
     997    float stepSize = params->data.F32[dim] * STEP_SIZE;
     998    //    float initialParam = params->data.F32[dim];
     999
     1000    if (stepSize == 0.0) {
     1001        stepSize = 1.0;
     1002    }
    9701003    a = b = c = params->data.F32[0];
    971     a = a * (1.0 - STEP_SIZE);
    972     c = c * (1.0 + STEP_SIZE);
    973 
    974     dummyParams->data.F32[0] = a;
    975     fa = func(dummyParams, coords);
    976     dummyParams->data.F32[0] = b;
    977     fb = func(dummyParams, coords);
    978     dummyParams->data.F32[0] = c;
    979     fc = func(dummyParams, coords);
     1004    a-= stepSize;
     1005    c+= stepSize;
     1006
     1007    params->data.F32[dim] = a;
     1008    fa = func(params, coords);
     1009    params->data.F32[dim] = b;
     1010    fb = func(params, coords);
     1011    params->data.F32[dim] = c;
     1012    fc = func(params, coords);
    9801013    if (fa < fb) {
    9811014        aDir = -1;
    … …  
    9971030            return(bracket);
    9981031        }
    999         a = a * (1.0 - STEP_SIZE);
    1000         c = c * (1.0 + STEP_SIZE);
    1001         dummyParams->data.F32[0] = a;
    1002         fa = func(dummyParams, coords);
    1003         dummyParams->data.F32[0] = c;
    1004         fc = func(dummyParams, coords);
     1032        stepSize*= (1.0 + stepSize);
     1033        a = a - stepSize;
     1034        c = c + stepSize;
     1035
     1036        params->data.F32[dim] = a;
     1037        fa = func(params, coords);
     1038        params->data.F32[dim] = c;
     1039        fc = func(params, coords);
     1040
     1041        //printf("HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    10051042
    10061043        if (fa < fb) {
    … …  
    10171054        if ((new_aDir == 1) && (aDir == -1)) {
    10181055            bracket->data.F32[0] = a;
    1019             bracket->data.F32[1] = bracket->data.F32[0] / (1.0 - STEP_SIZE);
    1020             bracket->data.F32[2] = bracket->data.F32[1] / (1.0 - STEP_SIZE);
     1056            bracket->data.F32[1] = b;
     1057            bracket->data.F32[2] = c;
    10211058            return(bracket);
    10221059        }
    10231060
    10241061        if ((new_cDir == 1) && (cDir == -1)) {
     1062            bracket->data.F32[0] = a;
     1063            bracket->data.F32[1] = b;
    10251064            bracket->data.F32[2] = c;
    1026             bracket->data.F32[1] = bracket->data.F32[0] / (1.0 + STEP_SIZE);
    1027             bracket->data.F32[0] = bracket->data.F32[1] / (1.0 + STEP_SIZE);
    10281065            return(bracket);
    10291066        }
    10301067        aDir = new_aDir;
    10311068        cDir = new_cDir;
    1032     }
    1033     psFree(dummyParams);
     1069        iter--;
     1070    }
    10341071    psFree(bracket);
    10351072    return(NULL);
    10361073}
    10371074
    1038 
    1039 
    10401075/******************************************************************************
    1041     This routine must minimize a function of exactly one parameter.
     1076    This routine takes as input an arbitrary function, and the parameter to
     1077    vary, and the line along which it must vary.  This function produces as
     1078    output a bracket [a, b, c] such that
     1079        f(param + b * line) is less than f(param + a * line) and
     1080                                         f(param + c * line).
     1081 
     1082    Algorithm: XXX completely ad hoc: start with the user-supplied starting
     1083    parameter and call that b.  Calculate a/c as a fractional amount
     1084    smaller/larger than b.  Repeat this process until a local minimum is
     1085    found.
     1086 *****************************************************************************/
     1087psVector *p_psDetermineBracket(psVector *params,
     1088                               psVector *line,
     1089                               const psVector *paramMask,
     1090                               const psArray *coords,
     1091                               psMinimizePowellFunc func)
     1092{
     1093    float a = 0.0;
     1094    float b = 0.0;
     1095    float c = 0.0;
     1096    float fa = 0.0;
     1097    float fb = 0.0;
     1098    float fc = 0.0;
     1099    int iter = 100;
     1100    float aDir = 0.0;
     1101    float cDir = 0.0;
     1102    float new_aDir = 0.0;
     1103    float new_cDir = 0.0;
     1104    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
     1105    float stepSize = STEP_SIZE;
     1106    psVector *tmp = NULL;
     1107    int i = 0;
     1108    int null = 0;
     1109
     1110    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1111            "---- p_psDetermineBracket() begin ----\n");
     1112
     1113    for (i=0;i<params->n;i++) {
     1114        if (paramMask->data.U8[i] == 0) {
     1115            if (line->data.F32[i] >= FLT_EPSILON) {
     1116                null = 1;
     1117            }
     1118        }
     1119    }
     1120    if (null == 0) {
     1121        psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
     1122                "p_psDetermineBracket() called with zero line vector.\n");
     1123        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1124                "---- p_psDetermineBracket() end (NULL) ----\n");
     1125        return(NULL);
     1126    }
     1127
     1128    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
     1129
     1130    b = 0;
     1131    a = -stepSize;
     1132    c = stepSize;
     1133
     1134    for (i=0;i<params->n;i++) {
     1135        if (paramMask->data.U8[i] == 0) {
     1136            tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1137        } else {
     1138            tmp->data.F32[i] = params->data.F32[i];
     1139        }
     1140    }
     1141    fa = func(tmp, coords);
     1142
     1143    for (i=0;i<params->n;i++) {
     1144        if (paramMask->data.U8[i] == 0) {
     1145            tmp->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
     1146        } else {
     1147            tmp->data.F32[i] = params->data.F32[i];
     1148        }
     1149    }
     1150    fb = func(tmp, coords);
     1151
     1152    for (i=0;i<params->n;i++) {
     1153        if (paramMask->data.U8[i] == 0) {
     1154            tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1155        } else {
     1156            tmp->data.F32[i] = params->data.F32[i];
     1157        }
     1158    }
     1159    fc = func(tmp, coords);
     1160
     1161    if (fa < fb) {
     1162        aDir = -1;
     1163    } else {
     1164        aDir = 1;
     1165    }
     1166
     1167    if (fc < fb) {
     1168        cDir = -1;
     1169    } else {
     1170        cDir = 1;
     1171    }
     1172
     1173    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1174            "HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1175    while (iter > 0) {
     1176        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1177                "psDetermineBracket(): iteration %d\n", iter);
     1178        if ((b < a) && (b < c)) {
     1179            bracket->data.F32[0] = a;
     1180            bracket->data.F32[1] = b;
     1181            bracket->data.F32[2] = c;
     1182            psFree(tmp);
     1183            psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1184                    "---- p_psDetermineBracket() end ----\n");
     1185            return(bracket);
     1186        }
     1187        stepSize*= (1.0 + stepSize);
     1188        a =- stepSize;
     1189        c =+ stepSize;
     1190
     1191        for (i=0;i<params->n;i++) {
     1192            if (paramMask->data.U8[i] == 0) {
     1193                tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1194            } else {
     1195                tmp->data.F32[i] = params->data.F32[i];
     1196            }
     1197        }
     1198        fa = func(tmp, coords);
     1199
     1200        for (i=0;i<params->n;i++) {
     1201            if (paramMask->data.U8[i] == 0) {
     1202                tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1203            } else {
     1204                tmp->data.F32[i] = params->data.F32[i];
     1205            }
     1206        }
     1207        fc = func(tmp, coords);
     1208
     1209        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1210                "HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1211
     1212        if (fa < fb) {
     1213            new_aDir = -1;
     1214        } else {
     1215            new_aDir = 1;
     1216        }
     1217
     1218        if (fc < fb) {
     1219            new_cDir = -1;
     1220        } else {
     1221            new_cDir = 1;
     1222        }
     1223        if ((new_aDir == 1) && (aDir == -1)) {
     1224            bracket->data.F32[0] = a;
     1225            bracket->data.F32[1] = b;
     1226            bracket->data.F32[2] = c;
     1227            psFree(tmp);
     1228            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1229                    "---- p_psDetermineBracket() end ----\n");
     1230            return(bracket);
     1231        }
     1232
     1233        if ((new_cDir == 1) && (cDir == -1)) {
     1234            bracket->data.F32[0] = a;
     1235            bracket->data.F32[1] = b;
     1236            bracket->data.F32[2] = c;
     1237            psFree(tmp);
     1238            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1239                    "---- p_psDetermineBracket() end ----\n");
     1240            return(bracket);
     1241        }
     1242        aDir = new_aDir;
     1243        cDir = new_cDir;
     1244        iter--;
     1245    }
     1246    psFree(tmp);
     1247    psFree(bracket);
     1248    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1249            "---- p_psDetermineBracket() end (NULL) ----\n");
     1250    return(NULL);
     1251}
     1252
     1253bool p_psDummyFunc(int dimNum,
     1254                   psVector *params,
     1255                   const psArray *coords,
     1256                   psMinimizePowellFunc func)
     1257{
     1258    return(true);
     1259}
     1260
     1261
     1262/******************************************************************************
     1263    This routine must minimize a possibly multi-dimensional function
     1264    (several parameters) along a single dimension.
    10421265 *****************************************************************************/
    10431266bool psMinimize1DFunc(psMinimization *min,
    10441267                      psVector *params,
     1268                      int dim,
    10451269                      const psArray *coords,
    10461270                      psMinimizePowellFunc func)
    10471271{
    10481272    psVector *bracket;
    1049     psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
    10501273    float a = 0.0;
    10511274    float b = 0.0;
    10521275    float c = 0.0;
    10531276    float n = 0.0;
     1277    float fa = 0.0;
    10541278    float fb = 0.0;
     1279    float fc = 0.0;
    10551280    float fn = 0.0;
    1056 
    1057     if (params->n != 1) {
    1058         psAbort(__func__, "params is larger than a single dimension.");
    1059     }
    1060     bracket = p_psDetermineBracket(params, coords, func);
    1061 
     1281    //    float initialParam = params->data.F32[dim];
     1282
     1283    bracket = p_psDetermineBracketOld(params, dim, coords, func);
     1284    if (bracket == NULL) {
     1285        psAbort(__func__, "Could not bracket minimum.");
     1286    }
    10621287
    10631288    min->iter = 0;
    1064     while (1) {
     1289    while (min->iter < min->maxIter) {
    10651290        min->iter++;
    1066 
     1291        //printf("psMinimize1DFunc(): iteration %d\n", min->iter);
    10671292        a = bracket->data.F32[0];
    10681293        b = bracket->data.F32[1];
    10691294        c = bracket->data.F32[2];
    1070         dummyParams->data.F32[0] = b;
    1071         fb = func(dummyParams, coords);
     1295
     1296        params->data.F32[dim] = a;
     1297        fa = func(params, coords);
     1298        params->data.F32[dim] = b;
     1299        fb = func(params, coords);
     1300        params->data.F32[dim] = c;
     1301        fc = func(params, coords);
     1302        //printf("Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
    10721303
    10731304        // We determine which is the biggest segment in [a,b,c] then split
    … …  
    10751306        if ((b-a) > (c-b)) {
    10761307            // This is the golden section formula
    1077             dummyParams->data.F32[0] = n = a + (0.69 * (b-a));
    1078             fn = func(dummyParams, coords);
     1308            params->data.F32[dim] = n = a + (0.69 * (b-a));
     1309            fn = func(params, coords);
    10791310            if (fn > fb) {
    10801311                // a = n, b = b, c = c
    … …  
    10861317            }
    10871318        } else {
    1088             dummyParams->data.F32[0] = b + (0.69 * (c-b));
    1089             fn = func(dummyParams, coords);
     1319            params->data.F32[dim] = n = b + (0.69 * (c-b));
     1320            fn = func(params, coords);
    10901321            if (fn > fb) {
    10911322                // a = a, b = b, c = n
    … …  
    10971328            }
    10981329        }
    1099     }
    1100 
    1101     psFree(dummyParams);
    1102 }
    1103 
     1330
     1331        if ((fabs(a-b) < min->tol) &&
     1332                (fabs(b-c) < min->tol)) {
     1333            //            psFree(bracket);
     1334            //  XXX: is this line correct?
     1335            params->data.F32[dim] = bracket->data.F32[1];
     1336            min->value = func(params, coords);
     1337            return(true);
     1338        }
     1339    }
     1340
     1341    //    psFree(bracket);
     1342    return(false);
     1343}
     1344
     1345/******************************************************************************
     1346    This routine must minimize a possibly multi-dimensional function
     1347    along a vector defined by line.
     1348 *****************************************************************************/
     1349float psLineMin(psMinimization *min,
     1350                psVector *params,
     1351                psVector *line,
     1352                const psVector *paramMask,
     1353                const psArray *coords,
     1354                psMinimizePowellFunc func)
     1355{
     1356    psVector *bracket;
     1357    float a = 0.0;
     1358    float b = 0.0;
     1359    float c = 0.0;
     1360    float n = 0.0;
     1361    float fa = 0.0;
     1362    float fb = 0.0;
     1363    float fc = 0.0;
     1364    float fn = 0.0;
     1365    float mul = 0.0;
     1366    psVector *tmpa = NULL;
     1367    psVector *tmpb = NULL;
     1368    psVector *tmpc = NULL;
     1369    psVector *tmpn = NULL;
     1370    int i = 0;
     1371    int null = 0;
     1372
     1373    psTrace(".psLib.dataManip.psLineMin", 4,
     1374            "---- psLineMin() begin ----\n");
     1375
     1376    for (i=0;i<params->n;i++) {
     1377        if (paramMask->data.U8[i] == 0) {
     1378            if (line->data.F32[i] >= FLT_EPSILON) {
     1379                null = 1;
     1380                psTrace(".psLib.dataManip.psLineMin", 4,
     1381                        "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
     1382            }
     1383        }
     1384    }
     1385    if (null == 0) {
     1386        psTrace(".psLib.dataManip.psLineMin", 2,
     1387                "psLineMin() called with zero line vector.\n");
     1388        return(0.0);
     1389    }
     1390
     1391    tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
     1392    tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
     1393    tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
     1394    tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
     1395
     1396
     1397    for (i=0;i<params->n;i++) {
     1398        psTrace(".psLib.dataManip.psLineMin", 6,
     1399                "params, paramMask, line [%d] is (%f %d %f)\n", i,
     1400                params->data.F32[i],
     1401                paramMask->data.U8[i],
     1402                line->data.F32[i]);
     1403    }
     1404
     1405    bracket = p_psDetermineBracket(params, line, paramMask, coords, func);
     1406    if (bracket == NULL) {
     1407        psAbort(__func__, "(1) Could not bracket minimum.");
     1408    }
     1409
     1410    min->iter = 0;
     1411    while (min->iter < min->maxIter) {
     1412        min->iter++;
     1413        psTrace(".psLib.dataManip.psLineMin", 6,
     1414                "psLineMin(): iteration %d\n", min->iter);
     1415
     1416        a = bracket->data.F32[0];
     1417        b = bracket->data.F32[1];
     1418        c = bracket->data.F32[2];
     1419
     1420        for (i=0;i<params->n;i++) {
     1421            if (paramMask->data.U8[i] == 0) {
     1422                tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1423                tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
     1424                tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1425            } else {
     1426                tmpa->data.F32[i] = params->data.F32[i];
     1427                tmpb->data.F32[i] = params->data.F32[i];
     1428                tmpc->data.F32[i] = params->data.F32[i];
     1429            }
     1430        }
     1431        fa = func(tmpa, coords);
     1432        fb = func(tmpb, coords);
     1433        fc = func(tmpc, coords);
     1434        psTrace(".psLib.dataManip.psLineMin", 6,
     1435                "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
     1436
     1437        // We determine which is the biggest segment in [a,b,c] then split
     1438        // that with the point n.
     1439        if ((b-a) > (c-b)) {
     1440            // This is the golden section formula
     1441            n = a + (0.69 * (b-a));
     1442            for (i=0;i<params->n;i++) {
     1443                tmpn->data.F32[i] = params->data.F32[i] + (n * line->data.F32[i]);
     1444            }
     1445            fn = func(tmpn, coords);
     1446
     1447            if (fn > fb) {
     1448                // a = n, b = b, c = c
     1449                bracket->data.F32[0] = n;
     1450            } else {
     1451                // a = a, b = n, c = b
     1452                bracket->data.F32[1] = n;
     1453                bracket->data.F32[2] = b;
     1454            }
     1455        } else {
     1456            n = b + (0.69 * (c-b));
     1457            for (i=0;i<params->n;i++) {
     1458                tmpn->data.F32[i] = params->data.F32[i] + (n * line->data.F32[i]);
     1459            }
     1460            fn = func(tmpn, coords);
     1461
     1462            if (fn > fb) {
     1463                // a = a, b = b, c = n
     1464                bracket->data.F32[2] = n;
     1465            } else {
     1466                // a = b, b = n, c = c
     1467                bracket->data.F32[0] = b;
     1468                bracket->data.F32[1] = n;
     1469            }
     1470        }
     1471        psTrace(".psLib.dataManip.psLineMin", 6,
     1472                "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
     1473
     1474        mul = bracket->data.F32[1];
     1475        if ((fabs(a-b) < min->tol) &&
     1476                (fabs(b-c) < min->tol)) {
     1477            for (i=0;i<params->n;i++) {
     1478                if (paramMask->data.U8[i] == 0) {
     1479                    params->data.F32[i]+= mul * line->data.F32[i];
     1480                }
     1481            }
     1482            min->value = func(params, coords);
     1483            psFree(bracket);
     1484            psFree(tmpa);
     1485            psFree(tmpb);
     1486            psFree(tmpc);
     1487            psFree(tmpn);
     1488            psTrace(".psLib.dataManip.psLineMin", 4,
     1489                    "---- psLineMin() end (%f) ----\n", mul);
     1490            return(mul);
     1491        }
     1492    }
     1493
     1494    psFree(bracket);
     1495    psFree(tmpa);
     1496    psFree(tmpb);
     1497    psFree(tmpc);
     1498    psFree(tmpn);
     1499
     1500    psTrace(".psLib.dataManip.psLineMin", 4,
     1501            "---- psLineMin() end (0.0) ----\n");
     1502    return(0.0);
     1503}
     1504
     1505
     1506/******************************************************************************
     1507    This routine must minimize a possibly multi-dimensional function.
     1508 *****************************************************************************/
    11041509bool psMinimizePowell(psMinimization *min,
    11051510                      psVector *params,
    … …  
    11081513                      psMinimizePowellFunc func)
    11091514{
    1110 
    1111 
    1112 
    1113 
    1114 
    1115 
    1116 
    1117 
    1118     return(true);
    1119 }
     1515    int numDims = params->n;
     1516    int i = 0;
     1517    int j = 0;
     1518    psVector **v = NULL;
     1519    psVector *pQP = NULL;
     1520    psVector *Q = NULL;
     1521    psVector *u = NULL;
     1522    psMinimization dummyMin;
     1523    float mul = 0.0;
     1524    float baseFuncVal = 0.0;
     1525    float currFuncVal = 0.0;
     1526    int biggestIter = 0;
     1527    float biggestDiff = 0.0;
     1528    float term1 = 0.0;
     1529    float term2 = 0.0;
     1530    float fqp = 0.0;
     1531    float diff = 0.0;
     1532
     1533    psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1534            "---- psMinimizePowell() begin ----\n");
     1535
     1536    pQP = psVectorAlloc(numDims, PS_TYPE_F32);
     1537    u = psVectorAlloc(numDims, PS_TYPE_F32);
     1538    Q = psVectorAlloc(numDims, PS_TYPE_F32);
     1539
     1540    // 1: Set v[i] to be the unit vectors for each dimension in params
     1541    v = (psVector **) psAlloc(params->n * sizeof(psVector **));
     1542    for (i=0;i<numDims;i++) {
     1543        v[i] = psVectorAlloc(numDims, PS_TYPE_F32);
     1544        for (j=0;j<numDims;j++) {
     1545            if (i == j) {
     1546                v[i]->data.F32[j] = 1.0;
     1547            } else {
     1548                v[i]->data.F32[j] = 0.0;
     1549            }
     1550        }
     1551    }
     1552
     1553    // 2: Set Q to be the initial params (P in the ADD)
     1554    for (i=0;i<numDims;i++) {
     1555        if (paramMask->data.U8[i] == 0) {
     1556            Q->data.F32[i] = params->data.F32[i];
     1557        }
     1558    }
     1559
     1560    min->iter = 0;
     1561    while (min->iter < min->maxIter) {
     1562        min->iter++;
     1563        psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1564                "psMinimizePowell() iteration %d\n", min->iter);
     1565
     1566        // 3: For each dimension in params, move Q only in the vector v[i] to
     1567        //    minimize the function.
     1568
     1569        baseFuncVal = func(Q, coords);
     1570        currFuncVal = baseFuncVal;
     1571        diff = 0.0;
     1572        biggestDiff = 0;
     1573        biggestIter = 0;
     1574        for (i=0;i<numDims;i++) {
     1575            if (paramMask->data.U8[i] == 0) {
     1576                mul = psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
     1577                if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
     1578                    biggestDiff = fabs(dummyMin.value - currFuncVal);
     1579                    biggestIter = i;
     1580                }
     1581                currFuncVal = dummyMin.value;
     1582            }
     1583        }
     1584
     1585        // 4: Set the vector u = Q - P
     1586        for (i=0;i<numDims;i++) {
     1587            if (paramMask->data.U8[i] == 0) {
     1588                u->data.F32[i] = Q->data.F32[i] - params->data.F32[i];
     1589
     1590                psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1591                        "WOW (%f %f %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]);
     1592
     1593            } else {
     1594                u->data.F32[i] = 0.0;
     1595            }
     1596        }
     1597
     1598        // 5: Move Q only in the direction u, and minimize the function.
     1599        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
     1600        psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1601                "u[0] is %f\n", u->data.F32[0]);
     1602
     1603        mul = psLineMin(min, params, u, paramMask, coords, func);
     1604
     1605        // 6:
     1606        if (min->value > currFuncVal) {
     1607            for (i=0;i<numDims;i++) {
     1608                psFree(v[i]);
     1609            }
     1610            psFree(v);
     1611            psFree(u);
     1612            psFree(pQP);
     1613            psFree(Q);
     1614            psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1615                    "---- psMinimizePowell() end (true) ----\n");
     1616            return(true);
     1617        }
     1618
     1619        for (i=0;i<numDims;i++) {
     1620            if (paramMask->data.U8[i] == 0) {
     1621                pQP->data.F32[i] = (2 * Q->data.F32[i]) - params->data.F32[i];
     1622            } else {
     1623                pQP->data.F32[i] = params->data.F32[i];
     1624            }
     1625        }
     1626        fqp = func(pQP, coords);
     1627
     1628        term1 = (baseFuncVal - currFuncVal) - biggestDiff;
     1629        term1*= term1;
     1630        term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp);
     1631        term2 = baseFuncVal - fqp;
     1632        term2*= term2 * biggestDiff;
     1633        if (term1 < term2) {
     1634            for (i=0;i<numDims;i++) {
     1635                if (paramMask->data.U8[i] == 0) {
     1636                    v[biggestIter]->data.F32[i] = u->data.F32[i];
     1637                }
     1638            }
     1639        }
     1640
     1641        // 7: Set P to Q
     1642        for (i=0;i<numDims;i++) {
     1643            if (paramMask->data.U8[i] == 0) {
     1644                params->data.F32[i] = Q->data.F32[i];
     1645            }
     1646        }
     1647
     1648        // 8: Go to step 3 until the change is less than some tolerance.
     1649        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
     1650            for (i=0;i<numDims;i++) {
     1651                psFree(v[i]);
     1652            }
     1653            psFree(v);
     1654            psFree(u);
     1655            psFree(pQP);
     1656            psFree(Q);
     1657            psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1658                    "---- psMinimizePowell() end (true) ----\n");
     1659            return(true);
     1660        }
     1661    }
     1662
     1663    for (i=0;i<numDims;i++) {
     1664        psFree(v[i]);
     1665    }
     1666    psFree(v);
     1667    psFree(u);
     1668    psFree(pQP);
     1669    psFree(Q);
     1670    psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1671            "---- psMinimizePowell() end (false) ----\n");
     1672    return(false);
     1673}

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

    r1629 r1719  
    105105typedef struct
    106106{
    107     const int maxIter;   ///< Convergence limit
    108     const float tol;   ///< Error Tolerance
     107    int maxIter;   ///< Convergence limit
     108    float tol;    ///< Error Tolerance
    109109    float value;   ///< Value of function at minimum
    110110    int iter;    ///< Number of iterations required
    … …  
    138138bool psMinimize1DFunc(psMinimization *min,
    139139                      psVector *params,
     140                      int dim,
    140141                      const psArray *coords,
    141142                      psMinimizePowellFunc func);
     143
     144float psLineMin(psMinimization *min,
     145                psVector *params,
     146                psVector *line,
     147                const psVector *paramMask,
     148                const psArray *coords,
     149                psMinimizePowellFunc func);
    142150
    143151

  • trunk/psLib/src/image/psImageStats.c

    r1448 r1719  
    1010*  @author George Gusciora, MHPCC
    1111*
    12 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-08-10 01:55:34 $
     12*  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-09-08 06:02:47 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    185185psPolynomial1D* *p_psCreateChebyshevPolys(int maxChebyPoly)
    186186{
    187     psPolynomial1D* *chebPolys = NULL;
     187    psPolynomial1D **chebPolys = NULL;
    188188    int i = 0;
    189189    int j = 0;
    190190
    191     chebPolys = (psPolynomial1D* *) psAlloc(maxChebyPoly * sizeof(psPolynomial1D* ));
     191    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
    192192    for (i = 0; i < maxChebyPoly; i++) {
    193         chebPolys[i] = psPolynomial1DAlloc(i + 1);
     193        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
    194194    }
    195195
    … …  
    333333    int i = 0;
    334334    int j = 0;
    335     float **sums = NULL;
     335    //    float **sums = NULL;
    336336    psPolynomial1D* *chebPolys = NULL;
    337337    int maxChebyPoly = 0;
    … …  
    340340    float polySum = 0.0;
    341341
    342     // Create the sums[][] data structure.  This
    343     // will hold the LHS of
    344     // equation
    345     // 29 in the ADD: sums[k][l] = SUM {
    346     // image(x,y) * Tk(x) * Tl(y) }
    347     sums = (float **)psAlloc(coeffs->nX * sizeof(float *));
    348     for (i = 0; i < coeffs->nX; i++) {
    349         sums[i] = (float *)psAlloc(coeffs->nY * sizeof(float));
    350     }
    351     for (i = 0; i < coeffs->nX; i++) {
    352         for (j = 0; j < coeffs->nY; j++) {
    353             sums[i][j] = 0.0;
    354         }
    355     }
     342    // Create the sums[][] data structure.  This will hold the LHS of
     343    // equation 29 in the ADD:
     344    //         sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
     345    //
     346    //    sums = (float **)psAlloc(coeffs->nX * sizeof(float *));
     347    //    for (i = 0; i < coeffs->nX; i++) {
     348    //        sums[i] = (float *)psAlloc(coeffs->nY * sizeof(float));
     349    //    }
     350    //    for (i = 0; i < coeffs->nX; i++) {
     351    //        for (j = 0; j < coeffs->nY; j++) {
     352    //            sums[i][j] = 0.0;
     353    //        }
     354    //    }
    356355
    357356    // We scale the pixel positions to values
    … …  
    375374                for (j = 0; j < coeffs->nY; j++) {
    376375                    polySum +=
    377                         psPolynomial1DEval
    378                         (rScalingFactors[x],
    379                          chebPolys[i]) *
    380                         psPolynomial1DEval(cScalingFactors[y], chebPolys[j]) * coeffs->coeff[i][j];
     376                        psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) *
     377                        psPolynomial1DEval(cScalingFactors[y], chebPolys[j]) *
     378                        coeffs->coeff[i][j];
    381379
    382380                }
    … …  
    394392
    395393    // Free some data
    396     for (i = 0; i < coeffs->nX; i++) {
    397         psFree(sums[i]);
    398     }
    399     psFree(sums);
     394    //    for (i = 0; i < coeffs->nX; i++) {
     395    //        psFree(sums[i]);
     396    //    }
     397    //    psFree(sums);
    400398    psFree(cScalingFactors);
    401399    psFree(rScalingFactors);

  • trunk/psLib/src/imageops/psImageStats.c

    r1448 r1719  
    1010*  @author George Gusciora, MHPCC
    1111*
    12 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-08-10 01:55:34 $
     12*  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-09-08 06:02:47 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    185185psPolynomial1D* *p_psCreateChebyshevPolys(int maxChebyPoly)
    186186{
    187     psPolynomial1D* *chebPolys = NULL;
     187    psPolynomial1D **chebPolys = NULL;
    188188    int i = 0;
    189189    int j = 0;
    190190
    191     chebPolys = (psPolynomial1D* *) psAlloc(maxChebyPoly * sizeof(psPolynomial1D* ));
     191    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
    192192    for (i = 0; i < maxChebyPoly; i++) {
    193         chebPolys[i] = psPolynomial1DAlloc(i + 1);
     193        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
    194194    }
    195195
    … …  
    333333    int i = 0;
    334334    int j = 0;
    335     float **sums = NULL;
     335    //    float **sums = NULL;
    336336    psPolynomial1D* *chebPolys = NULL;
    337337    int maxChebyPoly = 0;
    … …  
    340340    float polySum = 0.0;
    341341
    342     // Create the sums[][] data structure.  This
    343     // will hold the LHS of
    344     // equation
    345     // 29 in the ADD: sums[k][l] = SUM {
    346     // image(x,y) * Tk(x) * Tl(y) }
    347     sums = (float **)psAlloc(coeffs->nX * sizeof(float *));
    348     for (i = 0; i < coeffs->nX; i++) {
    349         sums[i] = (float *)psAlloc(coeffs->nY * sizeof(float));
    350     }
    351     for (i = 0; i < coeffs->nX; i++) {
    352         for (j = 0; j < coeffs->nY; j++) {
    353             sums[i][j] = 0.0;
    354         }
    355     }
     342    // Create the sums[][] data structure.  This will hold the LHS of
     343    // equation 29 in the ADD:
     344    //         sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
     345    //
     346    //    sums = (float **)psAlloc(coeffs->nX * sizeof(float *));
     347    //    for (i = 0; i < coeffs->nX; i++) {
     348    //        sums[i] = (float *)psAlloc(coeffs->nY * sizeof(float));
     349    //    }
     350    //    for (i = 0; i < coeffs->nX; i++) {
     351    //        for (j = 0; j < coeffs->nY; j++) {
     352    //            sums[i][j] = 0.0;
     353    //        }
     354    //    }
    356355
    357356    // We scale the pixel positions to values
    … …  
    375374                for (j = 0; j < coeffs->nY; j++) {
    376375                    polySum +=
    377                         psPolynomial1DEval
    378                         (rScalingFactors[x],
    379                          chebPolys[i]) *
    380                         psPolynomial1DEval(cScalingFactors[y], chebPolys[j]) * coeffs->coeff[i][j];
     376                        psPolynomial1DEval(rScalingFactors[x], chebPolys[i]) *
     377                        psPolynomial1DEval(cScalingFactors[y], chebPolys[j]) *
     378                        coeffs->coeff[i][j];
    381379
    382380                }
    … …  
    394392
    395393    // Free some data
    396     for (i = 0; i < coeffs->nX; i++) {
    397         psFree(sums[i]);
    398     }
    399     psFree(sums);
     394    //    for (i = 0; i < coeffs->nX; i++) {
     395    //        psFree(sums[i]);
     396    //    }
     397    //    psFree(sums);
    400398    psFree(cScalingFactors);
    401399    psFree(rScalingFactors);

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

    r1629 r1719  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-08-26 20:10:25 $
     11 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-08 06:02:23 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4343#include "psMinimize.h"
    4444#include "psMatrix.h"
     45#include "psTrace.h"
    4546
    4647/*****************************************************************************/
    … …  
    569570    gsl_vector *x;
    570571
     572    psTrace(".psLib.dataManip.psMinimize", 4,
     573            "---- psMinimize() begin ----\n");
     574
    571575    PS_CHECK_NULL_VECTOR(initialGuess);
    572576    PS_CHECK_EMPTY_VECTOR(initialGuess);
    … …  
    655659        }
    656660    }
     661    psTrace(".psLib.dataManip.psMinimize", 4,
     662            "---- psMinimize() end ----\n");
    657663    return (initialGuess);
    658664}
    … …  
    691697    float chiSqOld = 0.0;
    692698
     699    psTrace(".psLib.dataManip.psMinimizeChi2", 4,
     700            "---- psMinimizeChi2() begin ----\n");
     701
    693702    PS_CHECK_NULL_IMAGE(domain);
    694703    PS_CHECK_EMPTY_IMAGE(domain);
    … …  
    837846
    838847    // Bye bye.
     848    psTrace(".psLib.dataManip.psMinimizeChi2", 4,
     849            "---- psMinimizeChi2() end ----\n");
    839850    return (initialGuess);
    840851}
    … …  
    863874    psVector* xSums = NULL;
    864875
     876    psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
     877            "---- psVectorFitPolynomial1D() begin ----\n");
    865878    // printf("psVectorFitPolynomial1D()\n");
    866879    // for (i=0;i<x->n;i++) {
    … …  
    936949    psFree(xSums);
    937950
     951    psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
     952            "---- psVectorFitPolynomial1D() begin ----\n");
    938953    return (myPoly);
     954}
     955
     956psMinimization *psMinimizationAlloc(int maxIter,
     957                                    float tol)
     958{
     959    psMinimization *min = psAlloc(sizeof(psMinimization));
     960    min->maxIter = maxIter;
     961    min->tol = tol;
     962    min->value = 0.0;
     963    min->iter = 0;
     964    min->lastDelta = 0.0;
     965
     966    return(min);
    939967}
    940968
    941969#define STEP_SIZE 0.10
    942970/******************************************************************************
    943     This routine takes as input an arbitrary function of a single parameter.
    944     This function produces as output a bracket [a, b, c] such that f(b) is
    945     less than f(a) and f(b).
     971    This routine takes as input an arbitrary function, and the parameter to
     972    vary.  This function produces as output a bracket [a, b, c] such that
     973    f(b) is less than f(a) and f(b).
    946974 
    947975    Algorithm: XXX completely ad hoc: start with the user-supplied starting
    … …  
    950978    found.
    951979 *****************************************************************************/
    952 psVector *p_psDetermineBracket(psVector *params,
    953                                const psArray *coords,
    954                                psMinimizePowellFunc func)
     980psVector *p_psDetermineBracketOld(psVector *params,
     981                                  int dim,
     982                                  const psArray *coords,
     983                                  psMinimizePowellFunc func)
    955984{
    956985    float a = 0.0;
    … …  
    966995    float new_cDir = 0.0;
    967996    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    968     psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
    969 
     997    float stepSize = params->data.F32[dim] * STEP_SIZE;
     998    //    float initialParam = params->data.F32[dim];
     999
     1000    if (stepSize == 0.0) {
     1001        stepSize = 1.0;
     1002    }
    9701003    a = b = c = params->data.F32[0];
    971     a = a * (1.0 - STEP_SIZE);
    972     c = c * (1.0 + STEP_SIZE);
    973 
    974     dummyParams->data.F32[0] = a;
    975     fa = func(dummyParams, coords);
    976     dummyParams->data.F32[0] = b;
    977     fb = func(dummyParams, coords);
    978     dummyParams->data.F32[0] = c;
    979     fc = func(dummyParams, coords);
     1004    a-= stepSize;
     1005    c+= stepSize;
     1006
     1007    params->data.F32[dim] = a;
     1008    fa = func(params, coords);
     1009    params->data.F32[dim] = b;
     1010    fb = func(params, coords);
     1011    params->data.F32[dim] = c;
     1012    fc = func(params, coords);
    9801013    if (fa < fb) {
    9811014        aDir = -1;
    … …  
    9971030            return(bracket);
    9981031        }
    999         a = a * (1.0 - STEP_SIZE);
    1000         c = c * (1.0 + STEP_SIZE);
    1001         dummyParams->data.F32[0] = a;
    1002         fa = func(dummyParams, coords);
    1003         dummyParams->data.F32[0] = c;
    1004         fc = func(dummyParams, coords);
     1032        stepSize*= (1.0 + stepSize);
     1033        a = a - stepSize;
     1034        c = c + stepSize;
     1035
     1036        params->data.F32[dim] = a;
     1037        fa = func(params, coords);
     1038        params->data.F32[dim] = c;
     1039        fc = func(params, coords);
     1040
     1041        //printf("HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    10051042
    10061043        if (fa < fb) {
    … …  
    10171054        if ((new_aDir == 1) && (aDir == -1)) {
    10181055            bracket->data.F32[0] = a;
    1019             bracket->data.F32[1] = bracket->data.F32[0] / (1.0 - STEP_SIZE);
    1020             bracket->data.F32[2] = bracket->data.F32[1] / (1.0 - STEP_SIZE);
     1056            bracket->data.F32[1] = b;
     1057            bracket->data.F32[2] = c;
    10211058            return(bracket);
    10221059        }
    10231060
    10241061        if ((new_cDir == 1) && (cDir == -1)) {
     1062            bracket->data.F32[0] = a;
     1063            bracket->data.F32[1] = b;
    10251064            bracket->data.F32[2] = c;
    1026             bracket->data.F32[1] = bracket->data.F32[0] / (1.0 + STEP_SIZE);
    1027             bracket->data.F32[0] = bracket->data.F32[1] / (1.0 + STEP_SIZE);
    10281065            return(bracket);
    10291066        }
    10301067        aDir = new_aDir;
    10311068        cDir = new_cDir;
    1032     }
    1033     psFree(dummyParams);
     1069        iter--;
     1070    }
    10341071    psFree(bracket);
    10351072    return(NULL);
    10361073}
    10371074
    1038 
    1039 
    10401075/******************************************************************************
    1041     This routine must minimize a function of exactly one parameter.
     1076    This routine takes as input an arbitrary function, and the parameter to
     1077    vary, and the line along which it must vary.  This function produces as
     1078    output a bracket [a, b, c] such that
     1079        f(param + b * line) is less than f(param + a * line) and
     1080                                         f(param + c * line).
     1081 
     1082    Algorithm: XXX completely ad hoc: start with the user-supplied starting
     1083    parameter and call that b.  Calculate a/c as a fractional amount
     1084    smaller/larger than b.  Repeat this process until a local minimum is
     1085    found.
     1086 *****************************************************************************/
     1087psVector *p_psDetermineBracket(psVector *params,
     1088                               psVector *line,
     1089                               const psVector *paramMask,
     1090                               const psArray *coords,
     1091                               psMinimizePowellFunc func)
     1092{
     1093    float a = 0.0;
     1094    float b = 0.0;
     1095    float c = 0.0;
     1096    float fa = 0.0;
     1097    float fb = 0.0;
     1098    float fc = 0.0;
     1099    int iter = 100;
     1100    float aDir = 0.0;
     1101    float cDir = 0.0;
     1102    float new_aDir = 0.0;
     1103    float new_cDir = 0.0;
     1104    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
     1105    float stepSize = STEP_SIZE;
     1106    psVector *tmp = NULL;
     1107    int i = 0;
     1108    int null = 0;
     1109
     1110    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1111            "---- p_psDetermineBracket() begin ----\n");
     1112
     1113    for (i=0;i<params->n;i++) {
     1114        if (paramMask->data.U8[i] == 0) {
     1115            if (line->data.F32[i] >= FLT_EPSILON) {
     1116                null = 1;
     1117            }
     1118        }
     1119    }
     1120    if (null == 0) {
     1121        psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
     1122                "p_psDetermineBracket() called with zero line vector.\n");
     1123        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1124                "---- p_psDetermineBracket() end (NULL) ----\n");
     1125        return(NULL);
     1126    }
     1127
     1128    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
     1129
     1130    b = 0;
     1131    a = -stepSize;
     1132    c = stepSize;
     1133
     1134    for (i=0;i<params->n;i++) {
     1135        if (paramMask->data.U8[i] == 0) {
     1136            tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1137        } else {
     1138            tmp->data.F32[i] = params->data.F32[i];
     1139        }
     1140    }
     1141    fa = func(tmp, coords);
     1142
     1143    for (i=0;i<params->n;i++) {
     1144        if (paramMask->data.U8[i] == 0) {
     1145            tmp->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
     1146        } else {
     1147            tmp->data.F32[i] = params->data.F32[i];
     1148        }
     1149    }
     1150    fb = func(tmp, coords);
     1151
     1152    for (i=0;i<params->n;i++) {
     1153        if (paramMask->data.U8[i] == 0) {
     1154            tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1155        } else {
     1156            tmp->data.F32[i] = params->data.F32[i];
     1157        }
     1158    }
     1159    fc = func(tmp, coords);
     1160
     1161    if (fa < fb) {
     1162        aDir = -1;
     1163    } else {
     1164        aDir = 1;
     1165    }
     1166
     1167    if (fc < fb) {
     1168        cDir = -1;
     1169    } else {
     1170        cDir = 1;
     1171    }
     1172
     1173    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1174            "HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1175    while (iter > 0) {
     1176        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1177                "psDetermineBracket(): iteration %d\n", iter);
     1178        if ((b < a) && (b < c)) {
     1179            bracket->data.F32[0] = a;
     1180            bracket->data.F32[1] = b;
     1181            bracket->data.F32[2] = c;
     1182            psFree(tmp);
     1183            psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1184                    "---- p_psDetermineBracket() end ----\n");
     1185            return(bracket);
     1186        }
     1187        stepSize*= (1.0 + stepSize);
     1188        a =- stepSize;
     1189        c =+ stepSize;
     1190
     1191        for (i=0;i<params->n;i++) {
     1192            if (paramMask->data.U8[i] == 0) {
     1193                tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1194            } else {
     1195                tmp->data.F32[i] = params->data.F32[i];
     1196            }
     1197        }
     1198        fa = func(tmp, coords);
     1199
     1200        for (i=0;i<params->n;i++) {
     1201            if (paramMask->data.U8[i] == 0) {
     1202                tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1203            } else {
     1204                tmp->data.F32[i] = params->data.F32[i];
     1205            }
     1206        }
     1207        fc = func(tmp, coords);
     1208
     1209        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
     1210                "HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
     1211
     1212        if (fa < fb) {
     1213            new_aDir = -1;
     1214        } else {
     1215            new_aDir = 1;
     1216        }
     1217
     1218        if (fc < fb) {
     1219            new_cDir = -1;
     1220        } else {
     1221            new_cDir = 1;
     1222        }
     1223        if ((new_aDir == 1) && (aDir == -1)) {
     1224            bracket->data.F32[0] = a;
     1225            bracket->data.F32[1] = b;
     1226            bracket->data.F32[2] = c;
     1227            psFree(tmp);
     1228            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1229                    "---- p_psDetermineBracket() end ----\n");
     1230            return(bracket);
     1231        }
     1232
     1233        if ((new_cDir == 1) && (cDir == -1)) {
     1234            bracket->data.F32[0] = a;
     1235            bracket->data.F32[1] = b;
     1236            bracket->data.F32[2] = c;
     1237            psFree(tmp);
     1238            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1239                    "---- p_psDetermineBracket() end ----\n");
     1240            return(bracket);
     1241        }
     1242        aDir = new_aDir;
     1243        cDir = new_cDir;
     1244        iter--;
     1245    }
     1246    psFree(tmp);
     1247    psFree(bracket);
     1248    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     1249            "---- p_psDetermineBracket() end (NULL) ----\n");
     1250    return(NULL);
     1251}
     1252
     1253bool p_psDummyFunc(int dimNum,
     1254                   psVector *params,
     1255                   const psArray *coords,
     1256                   psMinimizePowellFunc func)
     1257{
     1258    return(true);
     1259}
     1260
     1261
     1262/******************************************************************************
     1263    This routine must minimize a possibly multi-dimensional function
     1264    (several parameters) along a single dimension.
    10421265 *****************************************************************************/
    10431266bool psMinimize1DFunc(psMinimization *min,
    10441267                      psVector *params,
     1268                      int dim,
    10451269                      const psArray *coords,
    10461270                      psMinimizePowellFunc func)
    10471271{
    10481272    psVector *bracket;
    1049     psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
    10501273    float a = 0.0;
    10511274    float b = 0.0;
    10521275    float c = 0.0;
    10531276    float n = 0.0;
     1277    float fa = 0.0;
    10541278    float fb = 0.0;
     1279    float fc = 0.0;
    10551280    float fn = 0.0;
    1056 
    1057     if (params->n != 1) {
    1058         psAbort(__func__, "params is larger than a single dimension.");
    1059     }
    1060     bracket = p_psDetermineBracket(params, coords, func);
    1061 
     1281    //    float initialParam = params->data.F32[dim];
     1282
     1283    bracket = p_psDetermineBracketOld(params, dim, coords, func);
     1284    if (bracket == NULL) {
     1285        psAbort(__func__, "Could not bracket minimum.");
     1286    }
    10621287
    10631288    min->iter = 0;
    1064     while (1) {
     1289    while (min->iter < min->maxIter) {
    10651290        min->iter++;
    1066 
     1291        //printf("psMinimize1DFunc(): iteration %d\n", min->iter);
    10671292        a = bracket->data.F32[0];
    10681293        b = bracket->data.F32[1];
    10691294        c = bracket->data.F32[2];
    1070         dummyParams->data.F32[0] = b;
    1071         fb = func(dummyParams, coords);
     1295
     1296        params->data.F32[dim] = a;
     1297        fa = func(params, coords);
     1298        params->data.F32[dim] = b;
     1299        fb = func(params, coords);
     1300        params->data.F32[dim] = c;
     1301        fc = func(params, coords);
     1302        //printf("Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
    10721303
    10731304        // We determine which is the biggest segment in [a,b,c] then split
    … …  
    10751306        if ((b-a) > (c-b)) {
    10761307            // This is the golden section formula
    1077             dummyParams->data.F32[0] = n = a + (0.69 * (b-a));
    1078             fn = func(dummyParams, coords);
     1308            params->data.F32[dim] = n = a + (0.69 * (b-a));
     1309            fn = func(params, coords);
    10791310            if (fn > fb) {
    10801311                // a = n, b = b, c = c
    … …  
    10861317            }
    10871318        } else {
    1088             dummyParams->data.F32[0] = b + (0.69 * (c-b));
    1089             fn = func(dummyParams, coords);
     1319            params->data.F32[dim] = n = b + (0.69 * (c-b));
     1320            fn = func(params, coords);
    10901321            if (fn > fb) {
    10911322                // a = a, b = b, c = n
    … …  
    10971328            }
    10981329        }
    1099     }
    1100 
    1101     psFree(dummyParams);
    1102 }
    1103 
     1330
     1331        if ((fabs(a-b) < min->tol) &&
     1332                (fabs(b-c) < min->tol)) {
     1333            //            psFree(bracket);
     1334            //  XXX: is this line correct?
     1335            params->data.F32[dim] = bracket->data.F32[1];
     1336            min->value = func(params, coords);
     1337            return(true);
     1338        }
     1339    }
     1340
     1341    //    psFree(bracket);
     1342    return(false);
     1343}
     1344
     1345/******************************************************************************
     1346    This routine must minimize a possibly multi-dimensional function
     1347    along a vector defined by line.
     1348 *****************************************************************************/
     1349float psLineMin(psMinimization *min,
     1350                psVector *params,
     1351                psVector *line,
     1352                const psVector *paramMask,
     1353                const psArray *coords,
     1354                psMinimizePowellFunc func)
     1355{
     1356    psVector *bracket;
     1357    float a = 0.0;
     1358    float b = 0.0;
     1359    float c = 0.0;
     1360    float n = 0.0;
     1361    float fa = 0.0;
     1362    float fb = 0.0;
     1363    float fc = 0.0;
     1364    float fn = 0.0;
     1365    float mul = 0.0;
     1366    psVector *tmpa = NULL;
     1367    psVector *tmpb = NULL;
     1368    psVector *tmpc = NULL;
     1369    psVector *tmpn = NULL;
     1370    int i = 0;
     1371    int null = 0;
     1372
     1373    psTrace(".psLib.dataManip.psLineMin", 4,
     1374            "---- psLineMin() begin ----\n");
     1375
     1376    for (i=0;i<params->n;i++) {
     1377        if (paramMask->data.U8[i] == 0) {
     1378            if (line->data.F32[i] >= FLT_EPSILON) {
     1379                null = 1;
     1380                psTrace(".psLib.dataManip.psLineMin", 4,
     1381                        "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
     1382            }
     1383        }
     1384    }
     1385    if (null == 0) {
     1386        psTrace(".psLib.dataManip.psLineMin", 2,
     1387                "psLineMin() called with zero line vector.\n");
     1388        return(0.0);
     1389    }
     1390
     1391    tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
     1392    tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
     1393    tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
     1394    tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
     1395
     1396
     1397    for (i=0;i<params->n;i++) {
     1398        psTrace(".psLib.dataManip.psLineMin", 6,
     1399                "params, paramMask, line [%d] is (%f %d %f)\n", i,
     1400                params->data.F32[i],
     1401                paramMask->data.U8[i],
     1402                line->data.F32[i]);
     1403    }
     1404
     1405    bracket = p_psDetermineBracket(params, line, paramMask, coords, func);
     1406    if (bracket == NULL) {
     1407        psAbort(__func__, "(1) Could not bracket minimum.");
     1408    }
     1409
     1410    min->iter = 0;
     1411    while (min->iter < min->maxIter) {
     1412        min->iter++;
     1413        psTrace(".psLib.dataManip.psLineMin", 6,
     1414                "psLineMin(): iteration %d\n", min->iter);
     1415
     1416        a = bracket->data.F32[0];
     1417        b = bracket->data.F32[1];
     1418        c = bracket->data.F32[2];
     1419
     1420        for (i=0;i<params->n;i++) {
     1421            if (paramMask->data.U8[i] == 0) {
     1422                tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
     1423                tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
     1424                tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
     1425            } else {
     1426                tmpa->data.F32[i] = params->data.F32[i];
     1427                tmpb->data.F32[i] = params->data.F32[i];
     1428                tmpc->data.F32[i] = params->data.F32[i];
     1429            }
     1430        }
     1431        fa = func(tmpa, coords);
     1432        fb = func(tmpb, coords);
     1433        fc = func(tmpc, coords);
     1434        psTrace(".psLib.dataManip.psLineMin", 6,
     1435                "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
     1436
     1437        // We determine which is the biggest segment in [a,b,c] then split
     1438        // that with the point n.
     1439        if ((b-a) > (c-b)) {
     1440            // This is the golden section formula
     1441            n = a + (0.69 * (b-a));
     1442            for (i=0;i<params->n;i++) {
     1443                tmpn->data.F32[i] = params->data.F32[i] + (n * line->data.F32[i]);
     1444            }
     1445            fn = func(tmpn, coords);
     1446
     1447            if (fn > fb) {
     1448                // a = n, b = b, c = c
     1449                bracket->data.F32[0] = n;
     1450            } else {
     1451                // a = a, b = n, c = b
     1452                bracket->data.F32[1] = n;
     1453                bracket->data.F32[2] = b;
     1454            }
     1455        } else {
     1456            n = b + (0.69 * (c-b));
     1457            for (i=0;i<params->n;i++) {
     1458                tmpn->data.F32[i] = params->data.F32[i] + (n * line->data.F32[i]);
     1459            }
     1460            fn = func(tmpn, coords);
     1461
     1462            if (fn > fb) {
     1463                // a = a, b = b, c = n
     1464                bracket->data.F32[2] = n;
     1465            } else {
     1466                // a = b, b = n, c = c
     1467                bracket->data.F32[0] = b;
     1468                bracket->data.F32[1] = n;
     1469            }
     1470        }
     1471        psTrace(".psLib.dataManip.psLineMin", 6,
     1472                "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
     1473
     1474        mul = bracket->data.F32[1];
     1475        if ((fabs(a-b) < min->tol) &&
     1476                (fabs(b-c) < min->tol)) {
     1477            for (i=0;i<params->n;i++) {
     1478                if (paramMask->data.U8[i] == 0) {
     1479                    params->data.F32[i]+= mul * line->data.F32[i];
     1480                }
     1481            }
     1482            min->value = func(params, coords);
     1483            psFree(bracket);
     1484            psFree(tmpa);
     1485            psFree(tmpb);
     1486            psFree(tmpc);
     1487            psFree(tmpn);
     1488            psTrace(".psLib.dataManip.psLineMin", 4,
     1489                    "---- psLineMin() end (%f) ----\n", mul);
     1490            return(mul);
     1491        }
     1492    }
     1493
     1494    psFree(bracket);
     1495    psFree(tmpa);
     1496    psFree(tmpb);
     1497    psFree(tmpc);
     1498    psFree(tmpn);
     1499
     1500    psTrace(".psLib.dataManip.psLineMin", 4,
     1501            "---- psLineMin() end (0.0) ----\n");
     1502    return(0.0);
     1503}
     1504
     1505
     1506/******************************************************************************
     1507    This routine must minimize a possibly multi-dimensional function.
     1508 *****************************************************************************/
    11041509bool psMinimizePowell(psMinimization *min,
    11051510                      psVector *params,
    … …  
    11081513                      psMinimizePowellFunc func)
    11091514{
    1110 
    1111 
    1112 
    1113 
    1114 
    1115 
    1116 
    1117 
    1118     return(true);
    1119 }
     1515    int numDims = params->n;
     1516    int i = 0;
     1517    int j = 0;
     1518    psVector **v = NULL;
     1519    psVector *pQP = NULL;
     1520    psVector *Q = NULL;
     1521    psVector *u = NULL;
     1522    psMinimization dummyMin;
     1523    float mul = 0.0;
     1524    float baseFuncVal = 0.0;
     1525    float currFuncVal = 0.0;
     1526    int biggestIter = 0;
     1527    float biggestDiff = 0.0;
     1528    float term1 = 0.0;
     1529    float term2 = 0.0;
     1530    float fqp = 0.0;
     1531    float diff = 0.0;
     1532
     1533    psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1534            "---- psMinimizePowell() begin ----\n");
     1535
     1536    pQP = psVectorAlloc(numDims, PS_TYPE_F32);
     1537    u = psVectorAlloc(numDims, PS_TYPE_F32);
     1538    Q = psVectorAlloc(numDims, PS_TYPE_F32);
     1539
     1540    // 1: Set v[i] to be the unit vectors for each dimension in params
     1541    v = (psVector **) psAlloc(params->n * sizeof(psVector **));
     1542    for (i=0;i<numDims;i++) {
     1543        v[i] = psVectorAlloc(numDims, PS_TYPE_F32);
     1544        for (j=0;j<numDims;j++) {
     1545            if (i == j) {
     1546                v[i]->data.F32[j] = 1.0;
     1547            } else {
     1548                v[i]->data.F32[j] = 0.0;
     1549            }
     1550        }
     1551    }
     1552
     1553    // 2: Set Q to be the initial params (P in the ADD)
     1554    for (i=0;i<numDims;i++) {
     1555        if (paramMask->data.U8[i] == 0) {
     1556            Q->data.F32[i] = params->data.F32[i];
     1557        }
     1558    }
     1559
     1560    min->iter = 0;
     1561    while (min->iter < min->maxIter) {
     1562        min->iter++;
     1563        psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1564                "psMinimizePowell() iteration %d\n", min->iter);
     1565
     1566        // 3: For each dimension in params, move Q only in the vector v[i] to
     1567        //    minimize the function.
     1568
     1569        baseFuncVal = func(Q, coords);
     1570        currFuncVal = baseFuncVal;
     1571        diff = 0.0;
     1572        biggestDiff = 0;
     1573        biggestIter = 0;
     1574        for (i=0;i<numDims;i++) {
     1575            if (paramMask->data.U8[i] == 0) {
     1576                mul = psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
     1577                if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
     1578                    biggestDiff = fabs(dummyMin.value - currFuncVal);
     1579                    biggestIter = i;
     1580                }
     1581                currFuncVal = dummyMin.value;
     1582            }
     1583        }
     1584
     1585        // 4: Set the vector u = Q - P
     1586        for (i=0;i<numDims;i++) {
     1587            if (paramMask->data.U8[i] == 0) {
     1588                u->data.F32[i] = Q->data.F32[i] - params->data.F32[i];
     1589
     1590                psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1591                        "WOW (%f %f %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]);
     1592
     1593            } else {
     1594                u->data.F32[i] = 0.0;
     1595            }
     1596        }
     1597
     1598        // 5: Move Q only in the direction u, and minimize the function.
     1599        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
     1600        psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1601                "u[0] is %f\n", u->data.F32[0]);
     1602
     1603        mul = psLineMin(min, params, u, paramMask, coords, func);
     1604
     1605        // 6:
     1606        if (min->value > currFuncVal) {
     1607            for (i=0;i<numDims;i++) {
     1608                psFree(v[i]);
     1609            }
     1610            psFree(v);
     1611            psFree(u);
     1612            psFree(pQP);
     1613            psFree(Q);
     1614            psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1615                    "---- psMinimizePowell() end (true) ----\n");
     1616            return(true);
     1617        }
     1618
     1619        for (i=0;i<numDims;i++) {
     1620            if (paramMask->data.U8[i] == 0) {
     1621                pQP->data.F32[i] = (2 * Q->data.F32[i]) - params->data.F32[i];
     1622            } else {
     1623                pQP->data.F32[i] = params->data.F32[i];
     1624            }
     1625        }
     1626        fqp = func(pQP, coords);
     1627
     1628        term1 = (baseFuncVal - currFuncVal) - biggestDiff;
     1629        term1*= term1;
     1630        term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp);
     1631        term2 = baseFuncVal - fqp;
     1632        term2*= term2 * biggestDiff;
     1633        if (term1 < term2) {
     1634            for (i=0;i<numDims;i++) {
     1635                if (paramMask->data.U8[i] == 0) {
     1636                    v[biggestIter]->data.F32[i] = u->data.F32[i];
     1637                }
     1638            }
     1639        }
     1640
     1641        // 7: Set P to Q
     1642        for (i=0;i<numDims;i++) {
     1643            if (paramMask->data.U8[i] == 0) {
     1644                params->data.F32[i] = Q->data.F32[i];
     1645            }
     1646        }
     1647
     1648        // 8: Go to step 3 until the change is less than some tolerance.
     1649        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
     1650            for (i=0;i<numDims;i++) {
     1651                psFree(v[i]);
     1652            }
     1653            psFree(v);
     1654            psFree(u);
     1655            psFree(pQP);
     1656            psFree(Q);
     1657            psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1658                    "---- psMinimizePowell() end (true) ----\n");
     1659            return(true);
     1660        }
     1661    }
     1662
     1663    for (i=0;i<numDims;i++) {
     1664        psFree(v[i]);
     1665    }
     1666    psFree(v);
     1667    psFree(u);
     1668    psFree(pQP);
     1669    psFree(Q);
     1670    psTrace(".psLib.dataManip.psMinimizePowell", 4,
     1671            "---- psMinimizePowell() end (false) ----\n");
     1672    return(false);
     1673}

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

    r1629 r1719  
    105105typedef struct
    106106{
    107     const int maxIter;   ///< Convergence limit
    108     const float tol;   ///< Error Tolerance
     107    int maxIter;   ///< Convergence limit
     108    float tol;    ///< Error Tolerance
    109109    float value;   ///< Value of function at minimum
    110110    int iter;    ///< Number of iterations required
    … …  
    138138bool psMinimize1DFunc(psMinimization *min,
    139139                      psVector *params,
     140                      int dim,
    140141                      const psArray *coords,
    141142                      psMinimizePowellFunc func);
     143
     144float psLineMin(psMinimization *min,
     145                psVector *params,
     146                psVector *line,
     147                const psVector *paramMask,
     148                const psArray *coords,
     149                psMinimizePowellFunc func);
    142150
    143151

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

    r1718 r1719  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 06:00:49 $
     9 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-08 06:01:21 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1313 */
    14 
    1514/*****************************************************************************/
    16 
    1715/*  INCLUDE FILES                                                            */
    18 
    1916/*****************************************************************************/
    2017#include <stdlib.h>
    … …  
    3734
    3835/*****************************************************************************/
    39 
    4036/* DEFINE STATEMENTS                                                         */
    41 
    4237/*****************************************************************************/
    4338
    … …  
    8378
    8479/*****************************************************************************/
    85 
    8680/* TYPE DEFINITIONS                                                          */
    87 
    8881/*****************************************************************************/
    8982static void polynomial1DFree(psPolynomial1D* myPoly);
    … …  
    9790
    9891/*****************************************************************************/
    99 
    10092/* GLOBAL VARIABLES                                                          */
    101 
    10293/*****************************************************************************/
    10394
    … …  
    10596
    10697/*****************************************************************************/
    107 
    10898/* FILE STATIC VARIABLES                                                     */
    109 
    11099/*****************************************************************************/
    111100
    … …  
    113102
    114103/*****************************************************************************/
    115 
    116104/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    117 
    118105/*****************************************************************************/
    119106
    … …  
    121108
    122109/*****************************************************************************/
    123 
    124110/*  FUNCTION IMPLEMENTATION - PUBLIC                                         */
    125 
    126111/*****************************************************************************/
    127112

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

    r1718 r1719  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 06:00:49 $
     9 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-08 06:01:21 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1313 */
    14 
    1514/*****************************************************************************/
    16 
    1715/*  INCLUDE FILES                                                            */
    18 
    1916/*****************************************************************************/
    2017#include <stdlib.h>
    … …  
    3734
    3835/*****************************************************************************/
    39 
    4036/* DEFINE STATEMENTS                                                         */
    41 
    4237/*****************************************************************************/
    4338
    … …  
    8378
    8479/*****************************************************************************/
    85 
    8680/* TYPE DEFINITIONS                                                          */
    87 
    8881/*****************************************************************************/
    8982static void polynomial1DFree(psPolynomial1D* myPoly);
    … …  
    9790
    9891/*****************************************************************************/
    99 
    10092/* GLOBAL VARIABLES                                                          */
    101 
    10293/*****************************************************************************/
    10394
    … …  
    10596
    10697/*****************************************************************************/
    107 
    10898/* FILE STATIC VARIABLES                                                     */
    109 
    11099/*****************************************************************************/
    111100
    … …  
    113102
    114103/*****************************************************************************/
    115 
    116104/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    117 
    118105/*****************************************************************************/
    119106
    … …  
    121108
    122109/*****************************************************************************/
    123 
    124110/*  FUNCTION IMPLEMENTATION - PUBLIC                                         */
    125 
    126111/*****************************************************************************/
    127112

  • trunk/psLib/src/sys/psTrace.h

    r1713 r1719  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-09-08 00:09:06 $
     12 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-09-08 06:03:26 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1919
    2020#define PS_UNKNOWN_TRACE_LEVEL -9999   // we don't know this name's level
    21 #define PS_DEFAULT_TRACE_LEVEL 0
     21#define PS_DEFAULT_TRACE_LEVEL -1
    2222
    2323/** \addtogroup LogTrace
Note: See TracChangeset for help on using the changeset viewer.