Changeset 7107 for trunk/psLib/src/math/psMinimizePolyFit.c
- Timestamp:
- May 10, 2006, 3:24:57 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psMinimizePolyFit.c (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psMinimizePolyFit.c
r7104 r7107 10 10 * @author EAM, IfA 11 11 * 12 * @version $Revision: 1.1 1$ $Name: not supported by cvs2svn $13 * @date $Date: 2006-05-10 1 1:38:55$12 * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2006-05-10 13:24:57 $ 14 14 * 15 15 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 1415 1415 } 1416 1416 1417 psImage *A = NULL; 1418 psVector *B = NULL; 1419 psF64 ***Sums = NULL; 1420 psF64 wt; 1421 psS32 nTerm; 1422 1423 psS32 nXterm = 1 + myPoly->nX; 1424 psS32 nYterm = 1 + myPoly->nY; 1425 psS32 nZterm = 1 + myPoly->nZ; 1426 nTerm = nXterm * nYterm * nZterm; 1427 1428 A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); 1429 B = psVectorAlloc(nTerm, PS_TYPE_F64); 1417 1418 int nXterm = 1 + myPoly->nX; // Number of x terms 1419 int nYterm = 1 + myPoly->nY; // Number of y terms 1420 int nZterm = 1 + myPoly->nZ; // Number of z terms 1421 int nTerm = nXterm * nYterm * nZterm; // Total number of terms 1422 int nData = x->n; // Number of data points 1423 psImage *A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); // Least-squares matrix 1424 psVector *B = psVectorAlloc(nTerm, PS_TYPE_F64); // Least-squares vector 1430 1425 B->n = B->nalloc; 1426 1431 1427 // Initialize data structures. 1432 1428 if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { … … 1439 1435 } 1440 1436 1441 // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ... 1437 // Dereference points for speed in the loop 1438 psF64 **matrix = A->data.F64; // Least-squares matrix 1439 psF64 *vector = B->data.F64; // Least-squares vector 1440 psF64 *xData = x->data.F64; // x 1441 psF64 *yData = y->data.F64; // y 1442 psF64 *zData = z->data.F64; // z 1443 psF64 *fData = f->data.F64; // f 1444 psF64 *fErrData = NULL; // Error in f 1445 if (fErr) { 1446 fErrData = fErr->data.F64; 1447 } 1448 psU8 *dataMask = NULL; // Mask for data 1449 if (mask) { 1450 dataMask = mask->data.U8; 1451 } 1452 psU8 ***termMask = myPoly->mask; // Mask for polynomial terms 1453 int nXYterm = nXterm * nYterm; // Multiplication of the numbers, to calculate the index 1442 1454 1443 1455 // Build the B and A data structs. 1444 for (psS32 k = 0; k < x->n; k++) { 1445 if ((mask != NULL) && (mask->data.U8[k] & maskValue)) { 1456 psF64 ***Sums = NULL; // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ... 1457 for (int k = 0; k < nData; k++) { 1458 if (dataMask && dataMask[k] & maskValue) { 1446 1459 continue; 1447 1460 } 1448 1461 1449 Sums = BuildSums3D(Sums, x->data.F64[k], y->data.F64[k], z->data.F64[k], nXterm, nYterm, nZterm); 1450 1462 Sums = BuildSums3D(Sums, xData[k], yData[k], zData[k], nXterm, nYterm, nZterm); 1463 1464 double wt; 1451 1465 if (fErr == NULL) { 1452 1466 wt = 1.0; 1453 1467 } else { 1454 1468 // this filters fErr == 0 values 1455 wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]); 1456 } 1457 1458 // we could skip half of the array and assign at the end 1459 // we must handle masked orders 1460 for (psS32 ix = 0; ix < nXterm; ix++) { 1461 for (psS32 iy = 0; iy < nYterm; iy++) { 1462 for (psS32 iz = 0; iz < nZterm; iz++) { 1463 if (myPoly->mask[ix][iy][iz]) { 1464 continue; 1465 } else { 1466 psS32 nx = ix + iy*nXterm + iz*nXterm*nYterm; 1467 B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt; 1468 } 1469 } 1470 } 1471 } 1472 1473 for (psS32 ix = 0; ix < nXterm; ix++) { 1474 for (psS32 iy = 0; iy < nYterm; iy++) { 1475 for (psS32 iz = 0; iz < nZterm; iz++) { 1476 if (myPoly->mask[ix][iy][iz]) 1477 continue; 1478 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm; 1479 for (psS32 jx = 0; jx < nXterm; jx++) { 1480 for (psS32 jy = 0; jy < nYterm; jy++) { 1481 for (psS32 jz = 0; jz < nZterm; jz++) { 1482 if (myPoly->mask[jx][jy][jz]) 1483 continue; 1484 psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm; 1485 A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt; 1486 } 1487 } 1488 } 1489 } 1490 } 1491 } 1492 } 1493 1494 for (psS32 ix = 0; ix < nXterm; ix++) { 1495 for (psS32 iy = 0; iy < nYterm; iy++) { 1496 for (psS32 iz = 0; iz < nZterm; iz++) { 1497 if (!myPoly->mask[ix][iy][iz]) 1469 wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErrData[k]); 1470 } 1471 1472 for (int i = 0; i < nTerm; i++) { 1473 int iz = i / nXYterm; // z index 1474 int iy = (i % nXYterm) / nXterm; // y index 1475 int ix = (i % nXYterm) % nXterm; // x index 1476 if (termMask[ix][iy][iz]) { 1477 continue; 1478 } 1479 1480 vector[i] += fData[k] * Sums[ix][iy][iz] * wt; 1481 matrix[i][i] += Sums[2*ix][2*iy][2*iz] * wt; 1482 for (int j = i + 1; j < nTerm; j++) { 1483 int jz = j / (nXterm * nYterm); // z index 1484 int jy = (j % (nXterm * nYterm)) / nXterm; // y index 1485 int jx = (j % (nXterm * nYterm)) % nXterm; // x index 1486 if (termMask[jx][jy][jz]) { 1498 1487 continue; 1499 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm; 1500 B->data.F64[nx] = 0; 1501 for (psS32 jx = 0; jx < nXterm; jx++) { 1502 for (psS32 jy = 0; jy < nYterm; jy++) { 1503 for (psS32 jz = 0; jz < nZterm; jz++) { 1504 psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm; 1505 A->data.F64[nx][ny] = (nx == ny) ? 1 : 0; 1506 } 1507 } 1508 } 1509 } 1510 } 1511 } 1488 } 1489 double value = Sums[ix+jx][iy+jy][iz+jz] * wt; 1490 matrix[i][j] += value; 1491 matrix[j][i] += value; 1492 } 1493 } 1494 } 1495 1496 // Free the sums 1497 for (psS32 ix = 0; ix < 2*nXterm; ix++) { 1498 for (psS32 iy = 0; iy < 2*nYterm; iy++) { 1499 psFree(Sums[ix][iy]); 1500 } 1501 psFree(Sums[ix]); 1502 } 1503 psFree(Sums); 1504 1512 1505 1513 1506 // XXX: rel10_ifa used psMatrixGJSolve(). However, this failed tests. So, I'm using psMatrixLUD(). … … 1518 1511 psFree(A); 1519 1512 psFree(B); 1520 1521 for (psS32 ix = 0; ix < 2*nXterm; ix++) {1522 for (psS32 iy = 0; iy < 2*nYterm; iy++) {1523 psFree(Sums[ix][iy]);1524 }1525 psFree(Sums[ix]);1526 }1527 psFree(Sums);1528 1513 psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n"); 1529 1514 return(NULL); … … 1578 1563 psFree(A); 1579 1564 psFree(B); 1580 1581 for (psS32 ix = 0; ix < 2*nXterm; ix++) {1582 for (psS32 iy = 0; iy < 2*nYterm; iy++) {1583 psFree(Sums[ix][iy]);1584 }1585 psFree(Sums[ix]);1586 }1587 psFree(Sums);1588 1565 1589 1566 psTrace(__func__, 4, "---- %s() end ----\n", __func__); … … 1961 1938 } 1962 1939 1963 // I think this is 1 dimension down 1964 psImage *A = NULL; 1965 psVector *B = NULL; 1966 psF64 ****Sums = NULL; 1967 psF64 wt; 1968 psS32 nTerm; 1969 1970 psS32 nXterm = 1 + myPoly->nX; 1971 psS32 nYterm = 1 + myPoly->nY; 1972 psS32 nZterm = 1 + myPoly->nZ; 1973 psS32 nTterm = 1 + myPoly->nZ; 1974 nTerm = nXterm * nYterm * nZterm * nTterm; 1975 1976 A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); 1977 B = psVectorAlloc(nTerm, PS_TYPE_F64); 1940 1941 int nXterm = 1 + myPoly->nX; // Number of x terms 1942 int nYterm = 1 + myPoly->nY; // Number of y terms 1943 int nZterm = 1 + myPoly->nZ; // Number of z terms 1944 int nTterm = 1 + myPoly->nT; // Number of t terms 1945 int nTerm = nXterm * nYterm * nZterm * nTterm; // Total number of terms 1946 int nData = x->n; // Number of data points 1947 psImage *A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); // Least-squares matrix 1948 psVector *B = psVectorAlloc(nTerm, PS_TYPE_F64); // Least-squares vector 1978 1949 B->n = B->nalloc; 1950 1979 1951 // Initialize data structures. 1980 1952 if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { … … 1987 1959 } 1988 1960 1989 // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ... 1961 // Dereference points for speed in the loop 1962 psF64 **matrix = A->data.F64; // Least-squares matrix 1963 psF64 *vector = B->data.F64; // Least-squares vector 1964 psF64 *xData = x->data.F64; // x 1965 psF64 *yData = y->data.F64; // y 1966 psF64 *zData = z->data.F64; // z 1967 psF64 *tData = t->data.F64; // t 1968 psF64 *fData = f->data.F64; // f 1969 psF64 *fErrData = NULL; // Error in f 1970 if (fErr) { 1971 fErrData = fErr->data.F64; 1972 } 1973 psU8 *dataMask = NULL; // Mask for data 1974 if (mask) { 1975 dataMask = mask->data.U8; 1976 } 1977 psU8 ****termMask = myPoly->mask; // Mask for polynomial terms 1978 int nXYZterm = nXterm * nYterm * nZterm; // Multiplication of the numbers, for calculating the index 1979 int nXYterm = nXterm * nYterm; // Multiplication of the numbers, for calculating the index 1990 1980 1991 1981 // Build the B and A data structs. 1992 for (psS32 k = 0; k < x->n; k++) { 1993 if ((mask != NULL) && (mask->data.U8[k] & maskValue)) { 1982 psF64 ****Sums = NULL; // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ... 1983 for (int k = 0; k < nData; k++) { 1984 if (dataMask && dataMask[k] & maskValue) { 1994 1985 continue; 1995 1986 } 1996 1987 1997 Sums = BuildSums4D(Sums, x->data.F64[k], y->data.F64[k], z->data.F64[k], t->data.F64[k], nXterm, nYterm, nZterm, nTterm); 1998 1988 Sums = BuildSums4D(Sums, xData[k], yData[k], zData[k], tData[k], nXterm, nYterm, nZterm, nTterm); 1989 1990 double wt; 1999 1991 if (fErr == NULL) { 2000 1992 wt = 1.0; 2001 1993 } else { 2002 1994 // this filters fErr == 0 values 2003 wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]); 2004 } 2005 2006 // we could skip half of the array and assign at the end 2007 // we must handle masked orders 2008 for (psS32 ix = 0; ix < nXterm; ix++) { 2009 for (psS32 iy = 0; iy < nYterm; iy++) { 2010 for (psS32 iz = 0; iz < nZterm; iz++) { 2011 for (psS32 it = 0; it < nTterm; it++) { 2012 if (myPoly->mask[ix][iy][iz][it]) 2013 continue; 2014 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; 2015 B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt; 2016 } 2017 } 2018 } 2019 } 2020 2021 for (psS32 ix = 0; ix < nXterm; ix++) { 2022 for (psS32 iy = 0; iy < nYterm; iy++) { 2023 for (psS32 iz = 0; iz < nZterm; iz++) { 2024 for (psS32 it = 0; it < nTterm; it++) { 2025 if (myPoly->mask[ix][iy][iz][it]) 2026 continue; 2027 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; 2028 for (psS32 jx = 0; jx < nXterm; jx++) { 2029 for (psS32 jy = 0; jy < nYterm; jy++) { 2030 for (psS32 jz = 0; jz < nZterm; jz++) { 2031 for (psS32 jt = 0; jt < nTterm; jt++) { 2032 if (myPoly->mask[jx][jy][jz][jt]) 2033 continue; 2034 psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm; 2035 A->data.F64[nx][ny]+= Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt; 2036 } 2037 } 2038 } 2039 } 2040 } 2041 } 2042 } 2043 } 2044 } 2045 2046 for (psS32 ix = 0; ix < nXterm; ix++) { 2047 for (psS32 iy = 0; iy < nYterm; iy++) { 2048 for (psS32 iz = 0; iz < nZterm; iz++) { 2049 for (psS32 it = 0; it < nTterm; it++) { 2050 if (!myPoly->mask[ix][iy][iz][it]) 2051 continue; 2052 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; 2053 B->data.F64[nx] = 0; 2054 for (psS32 jx = 0; jx < nXterm; jx++) { 2055 for (psS32 jy = 0; jy < nYterm; jy++) { 2056 for (psS32 jz = 0; jz < nZterm; jz++) { 2057 for (psS32 jt = 0; jt < nTterm; jt++) { 2058 psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm; 2059 A->data.F64[nx][ny] = (nx == ny) ? 1 : 0; 2060 } 2061 } 2062 } 2063 } 2064 } 2065 } 2066 } 2067 } 1995 wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErrData[k]); 1996 } 1997 1998 for (int i = 0; i < nTerm; i++) { 1999 int it = i / (nXYZterm); // t index 2000 int iz = (i % (nXYZterm)) / (nXYterm); // z index 2001 int iy = ((i % (nXYZterm)) % (nXYterm)) / nXterm; // y index 2002 int ix = ((i % (nXYZterm)) % (nXYterm)) % nXterm; // x index 2003 if (termMask[ix][iy][iz][it]) { 2004 continue; 2005 } 2006 2007 vector[i] += fData[k] * Sums[ix][iy][iz][it] * wt; 2008 matrix[i][i] += Sums[2*ix][2*iy][2*iz][2*it] * wt; 2009 for (int j = i + 1; j < nTerm; j++) { 2010 int jt = j / nXYZterm; // t index 2011 int jz = (j % nXYZterm) / nXYterm; // z index 2012 int jy = ((j % nXYZterm) % nXYterm) / nXterm; // y index 2013 int jx = ((j % nXYZterm) % nXYterm) % nXterm; // x index 2014 if (termMask[jx][jy][jz][jt]) { 2015 continue; 2016 } 2017 double value = Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt; 2018 matrix[i][j] += value; 2019 matrix[j][i] += value; 2020 } 2021 } 2022 } 2023 2024 // Free the sums 2025 for (int ix = 0; ix < 2*nXterm; ix++) { 2026 for (int iy = 0; iy < 2*nYterm; iy++) { 2027 for (int iz = 0; iz < 2*nZterm; iz++) { 2028 psFree(Sums[ix][iy][iz]); 2029 } 2030 psFree(Sums[ix][iy]); 2031 } 2032 psFree(Sums[ix]); 2033 } 2034 psFree(Sums); 2068 2035 2069 2036 … … 2072 2039 // does the solution in place 2073 2040 // The GaussJordan version was overflowing, so I'm using LUD. 2074 if ( psMatrixGJSolve(A, B)) {2041 if (!psMatrixGJSolve(A, B)) { 2075 2042 psFree(A); 2076 2043 psFree(B); 2077 for (psS32 ix = 0; ix < 2*nXterm; ix++) {2078 for (psS32 iy = 0; iy < 2*nYterm; iy++) {2079 for (psS32 iz = 0; iz < 2*nZterm; iz++) {2080 psFree(Sums[ix][iy][iz]);2081 }2082 psFree(Sums[ix][iy]);2083 }2084 psFree(Sums[ix]);2085 }2086 psFree(Sums);2087 2044 psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n"); 2088 2045 return(NULL); … … 2145 2102 psFree(B); 2146 2103 2147 for (psS32 ix = 0; ix < 2*nXterm; ix++) {2148 for (psS32 iy = 0; iy < 2*nYterm; iy++) {2149 for (psS32 iz = 0; iz < 2*nZterm; iz++) {2150 psFree(Sums[ix][iy][iz]);2151 }2152 psFree(Sums[ix][iy]);2153 }2154 psFree(Sums[ix]);2155 }2156 psFree(Sums);2157 2158 2104 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 2159 2105 return (myPoly);
Note:
See TracChangeset
for help on using the changeset viewer.
