IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1831


Ignore:
Timestamp:
Sep 17, 2004, 3:50:45 PM (22 years ago)
Author:
gusciora
Message:

Functions for psMinimizeLMChi2().

Location:
trunk/psLib/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r1823 r1831  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-17 02:14:52 $
     9 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-18 01:50:45 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    18091809p_psInterpolate1D(): This routine will take as input n-element floating
    18101810point arrays domain and range, and the x value, assumed to lie with the
    1811 domain vector.  It produces as output the n-order LaGrange interpolated
     1811domain vector.  It produces as output the (n-1)-order LaGrange interpolated
    18121812value of x.
    18131813 
    18141814XXX: do we error check for non-distinct domain values?
    18151815 *****************************************************************************/
    1816 float p_psFullInterpolate1DF32(float *domain,
     1816float p_ps1DFullInterpolateF32(float *domain,
    18171817                               float *range,
    18181818                               int n,
    … …  
    18221822    int m;
    18231823    static psVector *p = NULL;
    1824     psVectorRecycle(p, n, PS_TYPE_F32);
     1824    p = psVectorRecycle(p, n, PS_TYPE_F32);
    18251825    p_psMemSetPersistent(p, true);
     1826    p_psMemSetPersistent(p->data.F32, true);
     1827
     1828    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
     1829            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
     1830
     1831    for (i=0;i<n;i++) {
     1832        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1833                "domain/range is (%f %f)\n", domain[i], range[i]);
     1834    }
     1835
     1836    for (i=0;i<n;i++) {
     1837        p->data.F32[i] = range[i];
     1838        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1839                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
     1840
     1841    }
    18261842
    18271843    // From NR, during each iteration of the m loop, we are computing the
    18281844    // p_{i ... i+m} terms.
    1829     for (m=0;m<n;m++) {
     1845    for (m=1;m<n;m++) {
    18301846        for (i=0;i<n-m;i++) {
    1831             if (m == 0) {
    1832                 p->data.F32[i] = range[i];
    1833             } else {
    1834                 // From NR: we are computing P_{i ... i+m}
    1835                 p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
    1836                                   ((range[i]-x) * p->data.F32[i+1])) /
    1837                                  (domain[i] - domain[i+m]);
    1838             }
    1839         }
    1840     }
     1847            // From NR: we are computing P_{i ... i+m}
     1848            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
     1849                              ((domain[i]-x) * p->data.F32[i+1])) /
     1850                             (domain[i] - domain[i+m]);
     1851            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
     1852            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1853                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
     1854        }
     1855    }
     1856    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
     1857            "---- p_ps1DFullInterpolateF32() end ----\n");
     1858
    18411859    return(p->data.F32[0]);
    18421860}
    … …  
    18441862
    18451863// This is a
    1846 float p_psInterpolate1DF32(float *domain,
     1864float p_ps1DInterpolateF32(float *domain,
    18471865                           float *range,
    18481866                           int n,
    … …  
    18541872    int origin;
    18551873
     1874    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
     1875            "---- p_ps1DInterpolateF32() begin ----\n");
     1876
    18561877    binNum = VectorBinDisectF32(domain, n, x);
     1878
    18571879    if (0 == numIntPoints%2) {
    18581880        origin = binNum - ((numIntPoints/2) - 1);
    … …  
    18711893    }
    18721894
    1873     return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
     1895    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
     1896            "---- p_ps1DInterpolateF32() end ----\n");
     1897    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
    18741898}
    18751899
    … …  
    18821906XXX: This stuff does not work with a mask.
    18831907 *****************************************************************************/
    1884 float p_psInterpolate1D(psVector *domain,
    1885                         psVector *range,
    1886                         int order,
    1887                         psScalar *x)
    1888 {
     1908psScalar *p_psVectorInterpolate(psVector *domain,
     1909                                psVector *range,
     1910                                int order,
     1911                                psScalar *x)
     1912{
     1913    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1914            "---- p_psVectorInterpolate() begin ----\n");
     1915
    18891916    if (domain->type.type != range->type.type != x->type.type) {
    18901917        // XXX psError
    … …  
    18981925
    18991926    if (x->type.type == PS_TYPE_F32) {
    1900         return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
    1901                                     domain->n, order, x->data.F32));
     1927        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1928                "---- p_psVectorInterpolate() end ----\n");
     1929        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
     1930                             range->data.F32,
     1931                             domain->n,
     1932                             order,
     1933                             x->data.F32), PS_TYPE_F32));
    19021934    } else {
    19031935        // XXX psError: type not supported
    19041936    }
    19051937
    1906     return(-1.0);
     1938    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1939            "---- p_psVectorInterpolate() end ----\n");
     1940    return(NULL);
    19071941}
    19081942

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

    r1823 r1831  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-17 02:14:52 $
     14*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-18 01:50:45 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    409409                        psScalar *x);
    410410
     411psScalar *p_psVectorInterpolate(psVector *domain,
     412                                psVector *range,
     413                                int order,
     414                                psScalar *x);
     415
    411416/* \} */// End of MathGroup Functions
    412417

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

    r1819 r1831  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-16 19:30:24 $
     11 *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-18 01:50:45 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    844844    return (initialGuess);
    845845}
     846
    846847/******************************************************************************
    847848This routine will take an procedure which calculates an arbitrary function
    848849and it's derivative and minimize it.
    849850 
    850 GUS
    851  
    852851XXX: Do this:
    853852 After checking that all entries in the paramMask are 1 or 0, when
    … …  
    855854 
    856855     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    857  
    858  
    859856 *****************************************************************************/
    860857bool psMinimizeLM(psMinimization *min,
    … …  
    982979
    983980        oldValue = min->value;
    984         newValue = func(newDeriv, NRparams, coords);
     981        newValue = func(deriv, NRparams, coords);
    985982        psTrace(".psLib.dataManip.psMinimizeLM", 4,
    986983                "old/new values are (%f, %f)\n", oldValue, newValue);
    … …  
    10291026    psTrace(".psLib.dataManip.psMinimizeLM", 4,
    10301027            "---- psMinimizeLM() end (false) ----\n");
     1028    return(false);
     1029}
     1030
     1031/******************************************************************************
     1032This routine will take an procedure which calculates an arbitrary function
     1033and it's derivative and minimize the chi-squared match between that function
     1034at the specified coords and the specified value at those coords.
     1035 
     1036XXX: Do this:
     1037 After checking that all entries in the paramMask are 1 or 0, when
     1038 forming the A matrix from alpha, try this:
     1039 
     1040     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     1041 *****************************************************************************/
     1042bool psMinimizeLMChi2(psMinimization *min,
     1043                      psVector *params,
     1044                      const psVector *paramMask,
     1045                      psImage *covar,
     1046                      const psArray *coords,
     1047                      const psVector *value,
     1048                      psMinimizeLMChi2Func func)
     1049{
     1050    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1051            "---- psMinimizeLMChi2() begin ----\n");
     1052    int numData = value->n;
     1053    int numParams = params->n;
     1054    int i;
     1055    int j;
     1056    int k;
     1057    int l;
     1058    int n;
     1059    int p;
     1060    psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
     1061    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
     1062
     1063    psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
     1064    psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);
     1065    psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);
     1066
     1067    psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);
     1068    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1069    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1070
     1071    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
     1072    for (i=0;i<numData;i++) {
     1073        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
     1074    }
     1075
     1076    psVector *currValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
     1077    psVector *newValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
     1078
     1079    float currChi2 = 0.0;
     1080    float newChi2 = 0.0;
     1081    float lamda = 0.00005;
     1082
     1083    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1084            "min->maxIter is %d\n", min->maxIter);
     1085    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1086            "min->tol is %f\n", min->tol);
     1087
     1088    for (p=0;p<numParams;p++) {
     1089        origParams->data.F32[p] = params->data.F32[p];
     1090    }
     1091
     1092    min->lastDelta = HUGE;
     1093    min->lastDelta = 12345.0;
     1094    min->iter = 0;
     1095
     1096    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
     1097        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1098                "------------------------------------------------------\n");
     1099        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1100                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
     1101
     1102        //
     1103        // Calculate the current values and chi-squared of the function.
     1104        //
     1105        currChi2 = 0.0;
     1106        for (n=0;n<numData;n++) {
     1107            currValueVec->data.F32[n] = func(deriv[n],
     1108                                             params,
     1109                                             (psVector *) coords->data[n]);
     1110            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
     1111        }
     1112
     1113        for (p=0;p<numParams;p++) {
     1114            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1115                    "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
     1116        }
     1117        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1118                "Current chi-squared is (%f)\n", currChi2);
     1119
     1120        //
     1121        // Mask elements of the derivative for each data point.
     1122        //
     1123        for (p=0;p<numParams;p++) {
     1124            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
     1125                for (n=0;n<numData;n++) {
     1126                    (deriv[n])->data.F32[p] = 0.0;
     1127                }
     1128            }
     1129        }
     1130
     1131        //
     1132        // Calculate the BETA vector.
     1133        //
     1134        for (p=0;p<numParams;p++) {
     1135            beta->data.F64[p] = 0.0;
     1136            for (n=0;n<numData;n++) {
     1137                (beta->data.F64[p])+=
     1138                    (value->data.F32[n] - currValueVec->data.F32[n]) *
     1139                    (deriv[n])->data.F32[p];
     1140            }
     1141            // XXX: multiple by -1 here?
     1142            (beta->data.F64[p])*= -1.0;
     1143            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1144                    "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
     1145        }
     1146
     1147        //
     1148        // Calculate the ALPHA matrix.
     1149        //
     1150        for (k=0;k<numParams;k++) {
     1151            for (l=0;l<numParams;l++) {
     1152                alpha->data.F32[k][l] = 0.0;
     1153                for (n=0;n<numData;n++) {
     1154                    alpha->data.F32[k][l]+= (deriv[n])->data.F32[k] *
     1155                                            (deriv[n])->data.F32[l];
     1156                }
     1157            }
     1158        }
     1159
     1160        //
     1161        // Calculate the matrix A.
     1162        //
     1163        for (j=0;j<numParams;j++) {
     1164            for (k=0;k<numParams;k++) {
     1165                if (j == k) {
     1166                    A->data.F64[j][k] =
     1167                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
     1168                } else {
     1169                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
     1170                }
     1171            }
     1172        }
     1173        for (j=0;j<numParams;j++) {
     1174            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
     1175            for (k=0;k<numParams;k++) {
     1176                psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);
     1177            }
     1178            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
     1179        }
     1180
     1181        //
     1182        // Solve A * alpha = Beta
     1183        //
     1184        aOut = psMatrixLUD(aOut, perm, A);
     1185        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
     1186
     1187        //
     1188        // Mask any masked parameters.
     1189        //
     1190        for (i=0;i<numParams;i++) {
     1191            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1192                    "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
     1193            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
     1194                newParams->data.F32[i] = origParams->data.F32[i];
     1195            } else {
     1196                newParams->data.F32[i] = params->data.F32[i] -
     1197                                         (float) paramDeltasF64->data.F64[i];
     1198            }
     1199        }
     1200
     1201        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1202                "Calling func() with new parameters:\n");
     1203        for (i=0;i<numParams;i++) {
     1204            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1205                    "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
     1206        }
     1207
     1208
     1209        //
     1210        // Calculate new function values.
     1211        //
     1212        newChi2 = 0.0;
     1213        for (n=0;n<numData;n++) {
     1214            newValueVec->data.F32[n] = func(deriv[n], newParams,
     1215                                            (psVector *) coords->data[n]);
     1216            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
     1217        }
     1218        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1219                "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
     1220
     1221        //
     1222        // If the new chi-squared is lower, then keep it.
     1223        //
     1224        if (currChi2 > newChi2) {
     1225            min->lastDelta = currChi2 - newChi2;
     1226            min->value = newChi2;
     1227
     1228            // We already masked params.
     1229            for (i=0;i<numParams;i++) {
     1230                params->data.F32[i] = (float) newParams->data.F32[i];
     1231            }
     1232            lamda*= 0.1;
     1233        } else {
     1234            lamda*= 10.0;
     1235        }
     1236        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1237                "lamda is %f\n", lamda);
     1238        min->iter++;
     1239    }
     1240    psFree(beta);
     1241    psFree(perm);
     1242    psFree(paramDeltasF64);
     1243    psFree(origParams);
     1244    psFree(newParams);
     1245    psFree(alpha);
     1246    psFree(A);
     1247    psFree(aOut);
     1248    for (i=0;i<numData;i++) {
     1249        psFree(deriv[i]);
     1250    }
     1251    psFree(deriv);
     1252    psFree(currValueVec);
     1253    psFree(newValueVec);
     1254
     1255    if ((min->iter < min->maxIter) ||
     1256            (min->lastDelta <= min->tol)) {
     1257        return(true);
     1258    }
     1259
     1260    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1261            "---- psMinimizeLMChi2() end (false) ----\n");
    10311262    return(false);
    10321263}

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

    r1734 r1831  
    77 *  @author George Gusciora, MHPCC
    88 *
    9  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 23:32:03 $
     9 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-18 01:50:45 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    137137                                   const psArray *coords);
    138138
     139typedef float (*psMinimizeLMChi2Func) (psVector *deriv,
     140                                       const psVector *params,
     141                                       const psVector *coords);
     142
    139143bool psMinimizeLM(psMinimization *min,
    140144                  psImage *covar,
    … …  
    143147                  const psArray *coords,
    144148                  psMinimizeLMFunc func);
     149
     150bool psMinimizeLMChi2(psMinimization *min,
     151                      psVector *params,
     152                      const psVector *paramMask,
     153                      psImage *covar,
     154                      const psArray *coords,
     155                      const psVector *value,
     156                      psMinimizeLMChi2Func func);
    145157
    146158

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

    r1819 r1831  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-16 19:30:24 $
     11 *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-18 01:50:45 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    844844    return (initialGuess);
    845845}
     846
    846847/******************************************************************************
    847848This routine will take an procedure which calculates an arbitrary function
    848849and it's derivative and minimize it.
    849850 
    850 GUS
    851  
    852851XXX: Do this:
    853852 After checking that all entries in the paramMask are 1 or 0, when
    … …  
    855854 
    856855     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    857  
    858  
    859856 *****************************************************************************/
    860857bool psMinimizeLM(psMinimization *min,
    … …  
    982979
    983980        oldValue = min->value;
    984         newValue = func(newDeriv, NRparams, coords);
     981        newValue = func(deriv, NRparams, coords);
    985982        psTrace(".psLib.dataManip.psMinimizeLM", 4,
    986983                "old/new values are (%f, %f)\n", oldValue, newValue);
    … …  
    10291026    psTrace(".psLib.dataManip.psMinimizeLM", 4,
    10301027            "---- psMinimizeLM() end (false) ----\n");
     1028    return(false);
     1029}
     1030
     1031/******************************************************************************
     1032This routine will take an procedure which calculates an arbitrary function
     1033and it's derivative and minimize the chi-squared match between that function
     1034at the specified coords and the specified value at those coords.
     1035 
     1036XXX: Do this:
     1037 After checking that all entries in the paramMask are 1 or 0, when
     1038 forming the A matrix from alpha, try this:
     1039 
     1040     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     1041 *****************************************************************************/
     1042bool psMinimizeLMChi2(psMinimization *min,
     1043                      psVector *params,
     1044                      const psVector *paramMask,
     1045                      psImage *covar,
     1046                      const psArray *coords,
     1047                      const psVector *value,
     1048                      psMinimizeLMChi2Func func)
     1049{
     1050    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1051            "---- psMinimizeLMChi2() begin ----\n");
     1052    int numData = value->n;
     1053    int numParams = params->n;
     1054    int i;
     1055    int j;
     1056    int k;
     1057    int l;
     1058    int n;
     1059    int p;
     1060    psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
     1061    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
     1062
     1063    psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
     1064    psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);
     1065    psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);
     1066
     1067    psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);
     1068    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1069    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1070
     1071    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
     1072    for (i=0;i<numData;i++) {
     1073        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
     1074    }
     1075
     1076    psVector *currValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
     1077    psVector *newValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
     1078
     1079    float currChi2 = 0.0;
     1080    float newChi2 = 0.0;
     1081    float lamda = 0.00005;
     1082
     1083    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1084            "min->maxIter is %d\n", min->maxIter);
     1085    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1086            "min->tol is %f\n", min->tol);
     1087
     1088    for (p=0;p<numParams;p++) {
     1089        origParams->data.F32[p] = params->data.F32[p];
     1090    }
     1091
     1092    min->lastDelta = HUGE;
     1093    min->lastDelta = 12345.0;
     1094    min->iter = 0;
     1095
     1096    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
     1097        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1098                "------------------------------------------------------\n");
     1099        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1100                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
     1101
     1102        //
     1103        // Calculate the current values and chi-squared of the function.
     1104        //
     1105        currChi2 = 0.0;
     1106        for (n=0;n<numData;n++) {
     1107            currValueVec->data.F32[n] = func(deriv[n],
     1108                                             params,
     1109                                             (psVector *) coords->data[n]);
     1110            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
     1111        }
     1112
     1113        for (p=0;p<numParams;p++) {
     1114            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1115                    "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
     1116        }
     1117        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1118                "Current chi-squared is (%f)\n", currChi2);
     1119
     1120        //
     1121        // Mask elements of the derivative for each data point.
     1122        //
     1123        for (p=0;p<numParams;p++) {
     1124            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
     1125                for (n=0;n<numData;n++) {
     1126                    (deriv[n])->data.F32[p] = 0.0;
     1127                }
     1128            }
     1129        }
     1130
     1131        //
     1132        // Calculate the BETA vector.
     1133        //
     1134        for (p=0;p<numParams;p++) {
     1135            beta->data.F64[p] = 0.0;
     1136            for (n=0;n<numData;n++) {
     1137                (beta->data.F64[p])+=
     1138                    (value->data.F32[n] - currValueVec->data.F32[n]) *
     1139                    (deriv[n])->data.F32[p];
     1140            }
     1141            // XXX: multiple by -1 here?
     1142            (beta->data.F64[p])*= -1.0;
     1143            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1144                    "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
     1145        }
     1146
     1147        //
     1148        // Calculate the ALPHA matrix.
     1149        //
     1150        for (k=0;k<numParams;k++) {
     1151            for (l=0;l<numParams;l++) {
     1152                alpha->data.F32[k][l] = 0.0;
     1153                for (n=0;n<numData;n++) {
     1154                    alpha->data.F32[k][l]+= (deriv[n])->data.F32[k] *
     1155                                            (deriv[n])->data.F32[l];
     1156                }
     1157            }
     1158        }
     1159
     1160        //
     1161        // Calculate the matrix A.
     1162        //
     1163        for (j=0;j<numParams;j++) {
     1164            for (k=0;k<numParams;k++) {
     1165                if (j == k) {
     1166                    A->data.F64[j][k] =
     1167                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
     1168                } else {
     1169                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
     1170                }
     1171            }
     1172        }
     1173        for (j=0;j<numParams;j++) {
     1174            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
     1175            for (k=0;k<numParams;k++) {
     1176                psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);
     1177            }
     1178            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
     1179        }
     1180
     1181        //
     1182        // Solve A * alpha = Beta
     1183        //
     1184        aOut = psMatrixLUD(aOut, perm, A);
     1185        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
     1186
     1187        //
     1188        // Mask any masked parameters.
     1189        //
     1190        for (i=0;i<numParams;i++) {
     1191            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1192                    "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
     1193            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
     1194                newParams->data.F32[i] = origParams->data.F32[i];
     1195            } else {
     1196                newParams->data.F32[i] = params->data.F32[i] -
     1197                                         (float) paramDeltasF64->data.F64[i];
     1198            }
     1199        }
     1200
     1201        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1202                "Calling func() with new parameters:\n");
     1203        for (i=0;i<numParams;i++) {
     1204            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
     1205                    "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
     1206        }
     1207
     1208
     1209        //
     1210        // Calculate new function values.
     1211        //
     1212        newChi2 = 0.0;
     1213        for (n=0;n<numData;n++) {
     1214            newValueVec->data.F32[n] = func(deriv[n], newParams,
     1215                                            (psVector *) coords->data[n]);
     1216            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
     1217        }
     1218        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1219                "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
     1220
     1221        //
     1222        // If the new chi-squared is lower, then keep it.
     1223        //
     1224        if (currChi2 > newChi2) {
     1225            min->lastDelta = currChi2 - newChi2;
     1226            min->value = newChi2;
     1227
     1228            // We already masked params.
     1229            for (i=0;i<numParams;i++) {
     1230                params->data.F32[i] = (float) newParams->data.F32[i];
     1231            }
     1232            lamda*= 0.1;
     1233        } else {
     1234            lamda*= 10.0;
     1235        }
     1236        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1237                "lamda is %f\n", lamda);
     1238        min->iter++;
     1239    }
     1240    psFree(beta);
     1241    psFree(perm);
     1242    psFree(paramDeltasF64);
     1243    psFree(origParams);
     1244    psFree(newParams);
     1245    psFree(alpha);
     1246    psFree(A);
     1247    psFree(aOut);
     1248    for (i=0;i<numData;i++) {
     1249        psFree(deriv[i]);
     1250    }
     1251    psFree(deriv);
     1252    psFree(currValueVec);
     1253    psFree(newValueVec);
     1254
     1255    if ((min->iter < min->maxIter) ||
     1256            (min->lastDelta <= min->tol)) {
     1257        return(true);
     1258    }
     1259
     1260    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
     1261            "---- psMinimizeLMChi2() end (false) ----\n");
    10311262    return(false);
    10321263}

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

    r1734 r1831  
    77 *  @author George Gusciora, MHPCC
    88 *
    9  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-08 23:32:03 $
     9 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-18 01:50:45 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    137137                                   const psArray *coords);
    138138
     139typedef float (*psMinimizeLMChi2Func) (psVector *deriv,
     140                                       const psVector *params,
     141                                       const psVector *coords);
     142
    139143bool psMinimizeLM(psMinimization *min,
    140144                  psImage *covar,
    … …  
    143147                  const psArray *coords,
    144148                  psMinimizeLMFunc func);
     149
     150bool psMinimizeLMChi2(psMinimization *min,
     151                      psVector *params,
     152                      const psVector *paramMask,
     153                      psImage *covar,
     154                      const psArray *coords,
     155                      const psVector *value,
     156                      psMinimizeLMChi2Func func);
    145157
    146158

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

    r1823 r1831  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-17 02:14:52 $
     9 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-18 01:50:45 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    18091809p_psInterpolate1D(): This routine will take as input n-element floating
    18101810point arrays domain and range, and the x value, assumed to lie with the
    1811 domain vector.  It produces as output the n-order LaGrange interpolated
     1811domain vector.  It produces as output the (n-1)-order LaGrange interpolated
    18121812value of x.
    18131813 
    18141814XXX: do we error check for non-distinct domain values?
    18151815 *****************************************************************************/
    1816 float p_psFullInterpolate1DF32(float *domain,
     1816float p_ps1DFullInterpolateF32(float *domain,
    18171817                               float *range,
    18181818                               int n,
    … …  
    18221822    int m;
    18231823    static psVector *p = NULL;
    1824     psVectorRecycle(p, n, PS_TYPE_F32);
     1824    p = psVectorRecycle(p, n, PS_TYPE_F32);
    18251825    p_psMemSetPersistent(p, true);
     1826    p_psMemSetPersistent(p->data.F32, true);
     1827
     1828    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
     1829            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
     1830
     1831    for (i=0;i<n;i++) {
     1832        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1833                "domain/range is (%f %f)\n", domain[i], range[i]);
     1834    }
     1835
     1836    for (i=0;i<n;i++) {
     1837        p->data.F32[i] = range[i];
     1838        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1839                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
     1840
     1841    }
    18261842
    18271843    // From NR, during each iteration of the m loop, we are computing the
    18281844    // p_{i ... i+m} terms.
    1829     for (m=0;m<n;m++) {
     1845    for (m=1;m<n;m++) {
    18301846        for (i=0;i<n-m;i++) {
    1831             if (m == 0) {
    1832                 p->data.F32[i] = range[i];
    1833             } else {
    1834                 // From NR: we are computing P_{i ... i+m}
    1835                 p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
    1836                                   ((range[i]-x) * p->data.F32[i+1])) /
    1837                                  (domain[i] - domain[i+m]);
    1838             }
    1839         }
    1840     }
     1847            // From NR: we are computing P_{i ... i+m}
     1848            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
     1849                              ((domain[i]-x) * p->data.F32[i+1])) /
     1850                             (domain[i] - domain[i+m]);
     1851            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
     1852            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1853                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
     1854        }
     1855    }
     1856    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
     1857            "---- p_ps1DFullInterpolateF32() end ----\n");
     1858
    18411859    return(p->data.F32[0]);
    18421860}
    … …  
    18441862
    18451863// This is a
    1846 float p_psInterpolate1DF32(float *domain,
     1864float p_ps1DInterpolateF32(float *domain,
    18471865                           float *range,
    18481866                           int n,
    … …  
    18541872    int origin;
    18551873
     1874    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
     1875            "---- p_ps1DInterpolateF32() begin ----\n");
     1876
    18561877    binNum = VectorBinDisectF32(domain, n, x);
     1878
    18571879    if (0 == numIntPoints%2) {
    18581880        origin = binNum - ((numIntPoints/2) - 1);
    … …  
    18711893    }
    18721894
    1873     return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
     1895    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
     1896            "---- p_ps1DInterpolateF32() end ----\n");
     1897    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
    18741898}
    18751899
    … …  
    18821906XXX: This stuff does not work with a mask.
    18831907 *****************************************************************************/
    1884 float p_psInterpolate1D(psVector *domain,
    1885                         psVector *range,
    1886                         int order,
    1887                         psScalar *x)
    1888 {
     1908psScalar *p_psVectorInterpolate(psVector *domain,
     1909                                psVector *range,
     1910                                int order,
     1911                                psScalar *x)
     1912{
     1913    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1914            "---- p_psVectorInterpolate() begin ----\n");
     1915
    18891916    if (domain->type.type != range->type.type != x->type.type) {
    18901917        // XXX psError
    … …  
    18981925
    18991926    if (x->type.type == PS_TYPE_F32) {
    1900         return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
    1901                                     domain->n, order, x->data.F32));
     1927        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1928                "---- p_psVectorInterpolate() end ----\n");
     1929        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
     1930                             range->data.F32,
     1931                             domain->n,
     1932                             order,
     1933                             x->data.F32), PS_TYPE_F32));
    19021934    } else {
    19031935        // XXX psError: type not supported
    19041936    }
    19051937
    1906     return(-1.0);
     1938    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1939            "---- p_psVectorInterpolate() end ----\n");
     1940    return(NULL);
    19071941}
    19081942

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

    r1823 r1831  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-17 02:14:52 $
     14*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-18 01:50:45 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    409409                        psScalar *x);
    410410
     411psScalar *p_psVectorInterpolate(psVector *domain,
     412                                psVector *range,
     413                                int order,
     414                                psScalar *x);
     415
    411416/* \} */// End of MathGroup Functions
    412417

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

    r1823 r1831  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-17 02:14:52 $
     9 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-18 01:50:45 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    18091809p_psInterpolate1D(): This routine will take as input n-element floating
    18101810point arrays domain and range, and the x value, assumed to lie with the
    1811 domain vector.  It produces as output the n-order LaGrange interpolated
     1811domain vector.  It produces as output the (n-1)-order LaGrange interpolated
    18121812value of x.
    18131813 
    18141814XXX: do we error check for non-distinct domain values?
    18151815 *****************************************************************************/
    1816 float p_psFullInterpolate1DF32(float *domain,
     1816float p_ps1DFullInterpolateF32(float *domain,
    18171817                               float *range,
    18181818                               int n,
    … …  
    18221822    int m;
    18231823    static psVector *p = NULL;
    1824     psVectorRecycle(p, n, PS_TYPE_F32);
     1824    p = psVectorRecycle(p, n, PS_TYPE_F32);
    18251825    p_psMemSetPersistent(p, true);
     1826    p_psMemSetPersistent(p->data.F32, true);
     1827
     1828    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
     1829            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
     1830
     1831    for (i=0;i<n;i++) {
     1832        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1833                "domain/range is (%f %f)\n", domain[i], range[i]);
     1834    }
     1835
     1836    for (i=0;i<n;i++) {
     1837        p->data.F32[i] = range[i];
     1838        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1839                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
     1840
     1841    }
    18261842
    18271843    // From NR, during each iteration of the m loop, we are computing the
    18281844    // p_{i ... i+m} terms.
    1829     for (m=0;m<n;m++) {
     1845    for (m=1;m<n;m++) {
    18301846        for (i=0;i<n-m;i++) {
    1831             if (m == 0) {
    1832                 p->data.F32[i] = range[i];
    1833             } else {
    1834                 // From NR: we are computing P_{i ... i+m}
    1835                 p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
    1836                                   ((range[i]-x) * p->data.F32[i+1])) /
    1837                                  (domain[i] - domain[i+m]);
    1838             }
    1839         }
    1840     }
     1847            // From NR: we are computing P_{i ... i+m}
     1848            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
     1849                              ((domain[i]-x) * p->data.F32[i+1])) /
     1850                             (domain[i] - domain[i+m]);
     1851            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
     1852            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
     1853                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
     1854        }
     1855    }
     1856    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
     1857            "---- p_ps1DFullInterpolateF32() end ----\n");
     1858
    18411859    return(p->data.F32[0]);
    18421860}
    … …  
    18441862
    18451863// This is a
    1846 float p_psInterpolate1DF32(float *domain,
     1864float p_ps1DInterpolateF32(float *domain,
    18471865                           float *range,
    18481866                           int n,
    … …  
    18541872    int origin;
    18551873
     1874    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
     1875            "---- p_ps1DInterpolateF32() begin ----\n");
     1876
    18561877    binNum = VectorBinDisectF32(domain, n, x);
     1878
    18571879    if (0 == numIntPoints%2) {
    18581880        origin = binNum - ((numIntPoints/2) - 1);
    … …  
    18711893    }
    18721894
    1873     return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
     1895    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
     1896            "---- p_ps1DInterpolateF32() end ----\n");
     1897    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
    18741898}
    18751899
    … …  
    18821906XXX: This stuff does not work with a mask.
    18831907 *****************************************************************************/
    1884 float p_psInterpolate1D(psVector *domain,
    1885                         psVector *range,
    1886                         int order,
    1887                         psScalar *x)
    1888 {
     1908psScalar *p_psVectorInterpolate(psVector *domain,
     1909                                psVector *range,
     1910                                int order,
     1911                                psScalar *x)
     1912{
     1913    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1914            "---- p_psVectorInterpolate() begin ----\n");
     1915
    18891916    if (domain->type.type != range->type.type != x->type.type) {
    18901917        // XXX psError
    … …  
    18981925
    18991926    if (x->type.type == PS_TYPE_F32) {
    1900         return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
    1901                                     domain->n, order, x->data.F32));
     1927        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1928                "---- p_psVectorInterpolate() end ----\n");
     1929        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
     1930                             range->data.F32,
     1931                             domain->n,
     1932                             order,
     1933                             x->data.F32), PS_TYPE_F32));
    19021934    } else {
    19031935        // XXX psError: type not supported
    19041936    }
    19051937
    1906     return(-1.0);
     1938    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     1939            "---- p_psVectorInterpolate() end ----\n");
     1940    return(NULL);
    19071941}
    19081942

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

    r1823 r1831  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-17 02:14:52 $
     14*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-18 01:50:45 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    409409                        psScalar *x);
    410410
     411psScalar *p_psVectorInterpolate(psVector *domain,
     412                                psVector *range,
     413                                int order,
     414                                psScalar *x);
     415
    411416/* \} */// End of MathGroup Functions
    412417
Note: See TracChangeset for help on using the changeset viewer.