IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24086


Ignore:
Timestamp:
May 6, 2009, 2:19:32 PM (17 years ago)
Author:
eugene
Message:

explicitly call either GJ or LU solve functions; fixed errors on parameters (both solvers return inverse of A)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r21183 r24086  
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    1616 *
    17  *  XXX: psMatrixLUSolve() does not return error codes when the results are NANs.
     17 *  XXX: psMatrixLUSolution() does not return error codes when the results are NANs.
    1818 *
    1919 *  XXX: For clip-fit functions, what should we do if the mask is NULL?
    … …  
    403403    }
    404404
     405    bool status = false;
    405406    if (USE_GAUSS_JORDAN) {
    406         // GaussJordan version
    407         if (!psMatrixGJSolve(A, B)) {
    408             psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    409             goto escape_GJ;
    410         } else {
    411             // the first nTerm entries in B correspond directly to the desired
    412             // polynomial coefficients.  this is only true for the 1D case
    413             for (psS32 k = 0; k < numTerms; k++) {
    414                 myPoly->coeff[k] = B->data.F64[k];
    415                 myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
    416             }
    417             // The constant needs to be multiplied by 2, because it's half the a_0.
    418             myPoly->coeff[0] *= 2.0;
    419             myPoly->coeffErr[0] *= 2.0;
    420         }
     407        status = psMatrixGJSolve(A, B);
    421408    } else {
    422         // LUD version of the fit
    423         psImage *ALUD = NULL;
    424         psVector* outPerm = NULL;
    425         psVector* coeffs = NULL;
    426 
    427         ALUD = psImageAlloc(numTerms, numTerms, PS_TYPE_F64);
    428         ALUD = psMatrixLUD(ALUD, &outPerm, A);
    429         if (ALUD == NULL) {
    430             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
    431             goto escape_LUD;
    432         } else {
    433             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
    434             if (coeffs == NULL) {
    435                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
    436                 goto escape_LUD;
    437             } else {
    438                 for (psS32 k = 0; k < numTerms; k++) {
    439                     myPoly->coeff[k] = coeffs->data.F64[k];
    440                 }
    441             }
    442         }
    443         psFree(ALUD);
    444         psFree(coeffs);
    445         psFree(outPerm);
    446     }
     409        status = psMatrixLUSolve(A, B);
     410    }
     411    if (!status) {
     412        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
     413        goto escape;
     414    }
     415
     416    // the first nTerm entries in B correspond directly to the desired
     417    // polynomial coefficients.  this is only true for the 1D case
     418    for (psS32 k = 0; k < numTerms; k++) {
     419        myPoly->coeff[k] = B->data.F64[k];
     420        myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
     421    }
     422    // The constant needs to be multiplied by 2, because it's half the a_0.
     423    myPoly->coeff[0] *= 2.0;
     424    myPoly->coeffErr[0] *= 2.0;
     425
    447426    psFree(A);
    448427    psFree(B);
    449428    return true;
    450429
    451 escape_LUD:
    452     // XXX drop the LUD version!
    453     psFree(A);
    454     psFree(B);
    455     return false;
    456 
    457 escape_GJ:
     430escape:
    458431    psFree(A);
    459432    psFree(B);
    … …  
    628601    }
    629602
    630     // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
     603    bool status = false;
    631604    if (USE_GAUSS_JORDAN) {
    632         // GaussJordan version
    633         if (!psMatrixGJSolve(A, B)) {
    634             psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    635             goto escape_GJ;
    636         } else {
    637             // the first nTerm entries in B correspond directly to the desired
    638             // polynomial coefficients.  this is only true for the 1D case
    639             for (psS32 k = 0; k < nTerm; k++) {
    640                 if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
    641                 myPoly->coeff[k] = B->data.F64[k];
    642                 myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
    643             }
    644         }
     605        status = psMatrixGJSolve(A, B);
    645606    } else {
    646         // LUD version of the fit
    647         psImage *ALUD = NULL;
    648         psVector* outPerm = NULL;
    649         psVector* coeffs = NULL;
    650 
    651         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    652         ALUD = psMatrixLUD(ALUD, &outPerm, A);
    653         if (ALUD == NULL) {
    654             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
    655             goto escape_LUD;
    656         } else {
    657             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
    658             if (coeffs == NULL) {
    659                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
    660                 goto escape_LUD;
    661             } else {
    662                 for (psS32 k = 0; k < nTerm; k++) {
    663                     if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
    664                     myPoly->coeff[k] = coeffs->data.F64[k];
    665                     // XXX LUD does not give inverse of A
    666                 }
    667             }
    668         }
    669         psFree(ALUD);
    670         psFree(coeffs);
    671         psFree(outPerm);
    672     }
     607        status = psMatrixLUSolve(A, B);
     608    }
     609    if (!status) {
     610        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
     611        goto escape;
     612    }
     613
     614    // the first nTerm entries in B correspond directly to the desired
     615    // polynomial coefficients.  this is only true for the 1D case
     616    for (psS32 k = 0; k < nTerm; k++) {
     617        if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
     618        myPoly->coeff[k] = B->data.F64[k];
     619        myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
     620    }
     621
    673622    psFree(A);
    674623    psFree(B);
    675 
    676     psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
    677624    return true;
    678625
    679 escape_LUD:
    680     psFree(A);
    681     psFree(B);
    682     return false;
    683 
    684 escape_GJ:
     626escape:
    685627    psFree(A);
    686628    psFree(B);
    687629    return false;
    688630}
    689 
    690631
    691632/******************************************************************************
    … …  
    11291070    }
    11301071
    1131     if (!psMatrixGJSolve(A, B)) {
    1132         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    1133         psFree(A);
    1134         psFree(B);
    1135         return false;
    1136     }
    1137 
    1138     // select the appropriate solution entries
     1072    bool status = false;
     1073    if (USE_GAUSS_JORDAN) {
     1074        status = psMatrixGJSolve(A, B);
     1075    } else {
     1076        status = psMatrixLUSolve(A, B);
     1077    }
     1078    if (!status) {
     1079        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
     1080        goto escape;
     1081    }
     1082
     1083    // select the appropriate solution entries (retain the incoming values if masked on the fit)
    11391084    for (int i = 0; i < nTerm; i++) {
    1140         int l = i / nYterm;         // x index
    1141         int m = i % nYterm;         // y index
    1142 
    1143         // retain the incoming values if masked on the fit
    1144         if (coeffMask[l][m] & PS_POLY_MASK_FIT) continue;
    1145         myPoly->coeff[l][m] = B->data.F64[i];
    1146         myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
     1085        int ix = i / nYterm;         // x index
     1086        int iy = i % nYterm;         // y index
     1087        if (coeffMask[ix][iy] & PS_POLY_MASK_FIT) continue;
     1088        myPoly->coeff[ix][iy] = B->data.F64[i];
     1089        myPoly->coeffErr[ix][iy] = sqrt(A->data.F64[i][i]);
    11471090    }
    11481091    psFree(A);
    11491092    psFree(B);
    1150 
    1151     psTrace("psLib.math", 6, "---- %s() end ----\n", __func__);
    11521093    return true;
     1094
     1095escape:
     1096    psFree (A);
     1097    psFree (B);
     1098    return false;
    11531099}
    11541100
    … …  
    15231469
    15241470
    1525     // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
     1471    bool status = false;
    15261472    if (USE_GAUSS_JORDAN) {
    1527         // does the solution in place
    1528         // The matrices were overflowing, so I switched to LUD.
    1529         if (!psMatrixGJSolve(A, B)) {
    1530             psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
    1531             goto escape_GJ;
    1532         }
    1533         // select the appropriate solution entries
    1534         for (int i = 0; i < nTerm; i++) {
    1535             int ix = i / nYZterm; // x index
    1536             int iy = (i % nYZterm) / nZterm; // y index
    1537             int iz = (i % nYZterm) % nZterm; // z index
    1538             if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
    1539             myPoly->coeff[ix][iy][iz] = B->data.F64[i];
    1540             myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
    1541         }
     1473        status = psMatrixGJSolve(A, B);
    15421474    } else {
    1543         // LUD version of the fit
    1544         psImage *ALUD = NULL;
    1545         psVector* outPerm = NULL;
    1546         psVector* coeffs = NULL;
    1547 
    1548         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    1549         ALUD = psMatrixLUD(ALUD, &outPerm, A);
    1550         if (ALUD == NULL) {
    1551             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
    1552             goto escape_LUD;
    1553         } else {
    1554             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
    1555             if (coeffs == NULL) {
    1556                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
    1557                 goto escape_LUD;
    1558             } else {
    1559                 // select the appropriate solution entries
    1560                 for (psS32 ix = 0; ix < nXterm; ix++) {
    1561                     for (psS32 iy = 0; iy < nYterm; iy++) {
    1562                         for (psS32 iz = 0; iz < nZterm; iz++) {
    1563                             psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
    1564                             if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
    1565                             myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
    1566                             // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
    1567                             // XXX this is wrong: A is not replaced with inverse(A), as for GaussJordan
    1568                         }
    1569                     }
    1570                 }
    1571             }
    1572         }
    1573 
    1574         psFree(ALUD);
    1575         psFree(coeffs);
    1576         psFree(outPerm);
     1475        status = psMatrixLUSolve(A, B);
     1476    }
     1477    if (!status) {
     1478        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
     1479        goto escape;
     1480    }
     1481
     1482    // select the appropriate solution entries
     1483    for (int i = 0; i < nTerm; i++) {
     1484        int ix = i / nYZterm; // x index
     1485        int iy = (i % nYZterm) / nZterm; // y index
     1486        int iz = (i % nYZterm) % nZterm; // z index
     1487        if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
     1488        myPoly->coeff[ix][iy][iz] = B->data.F64[i];
     1489        myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
    15771490    }
    15781491    psFree(A);
    15791492    psFree(B);
    1580 
    1581     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
    15821493    return true;
    15831494
    1584 escape_LUD:
    1585     psFree(A);
    1586     psFree(B);
    1587     return false;
    1588 
    1589 escape_GJ:
    1590 
     1495escape:
    15911496    psFree(A);
    15921497    psFree(B);
    … …  
    19861891    }
    19871892
    1988     // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
     1893    bool status = false;
    19891894    if (USE_GAUSS_JORDAN) {
    1990         // does the solution in place
    1991         // The GaussJordan version was overflowing, so I'm using LUD.
    1992         if (!psMatrixGJSolve(A, B)) {
    1993             psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
    1994             goto escape_GJ;
    1995         }
    1996 
    1997         // select the appropriate solution entries
    1998         for (int i = 0; i < nTerm; i++) {
    1999             int ix = i / nYZTterm; // x index
    2000             int iy = (i % nYZTterm) / nZTterm; // y index
    2001             int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
    2002             int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
    2003             if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
    2004             myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
    2005             myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
    2006         }
     1895        status = psMatrixGJSolve(A, B);
    20071896    } else {
    2008         // LUD version of the fit
    2009         psImage *ALUD = NULL;
    2010         psVector* outPerm = NULL;
    2011         psVector* coeffs = NULL;
    2012 
    2013         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    2014         ALUD = psMatrixLUD(ALUD, &outPerm, A);
    2015         if (ALUD == NULL) {
    2016             psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
    2017             goto escape_LUD;
    2018         } else {
    2019             coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
    2020             if (coeffs == NULL) {
    2021                 psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
    2022                 goto escape_LUD;
    2023             } else {
    2024                 // select the appropriate solution entries
    2025                 for (psS32 ix = 0; ix < nXterm; ix++) {
    2026                     for (psS32 iy = 0; iy < nYterm; iy++) {
    2027                         for (psS32 iz = 0; iz < nZterm; iz++) {
    2028                             for (psS32 it = 0; it < nTterm; it++) {
    2029                                 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    2030                                 if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
    2031                                 myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
    2032                                 // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
    2033                                 // XXX this is wrong: LUD does not supply inverse(A)
    2034                             }
    2035                         }
    2036                     }
    2037                 }
    2038             }
    2039         }
    2040         psFree(ALUD);
    2041         psFree(coeffs);
    2042         psFree(outPerm);
     1897        status = psMatrixLUSolve(A, B);
     1898    }
     1899    if (!status) {
     1900        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
     1901        goto escape;
     1902    }
     1903
     1904    // select the appropriate solution entries
     1905    for (int i = 0; i < nTerm; i++) {
     1906        int ix = i / nYZTterm; // x index
     1907        int iy = (i % nYZTterm) / nZTterm; // y index
     1908        int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
     1909        int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
     1910        if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
     1911        myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
     1912        myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
    20431913    }
    20441914    psFree(A);
    20451915    psFree(B);
    2046 
    2047     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
    20481916    return true;
    20491917
    2050 escape_LUD:
    2051     psFree(A);
    2052     psFree(B);
    2053     return false;
    2054 
    2055 escape_GJ:
     1918escape:
    20561919    psFree(A);
    20571920    psFree(B);
Note: See TracChangeset for help on using the changeset viewer.