Changeset 24087

Timestamp:

May 6, 2009, 2:19:50 PM (17 years ago)

Author:

eugene

Message:

add a test for an example ill-conditioned matrix

Location:

trunk/psLib/test/math

Files:

• Makefile.am (modified) (1 diff)
• data (added)
• data/Agj.fits (added)
• data/Bgj.fits (added)
• data/input (added)
• tap_psMatrix08.c (modified) (3 diffs)

trunk/psLib/test/math/Makefile.am

@@ -48 +48 @@
         tap_psMatrixVectorArithmetic01 \
         tap_psMatrixVectorArithmetic04 \
-        tap_psRandom \
         tap_psMinimizePowell \
         tap_psSpline1D \
         tap_psPolynomialMD
+
+#       tap_psRandom
 
 if BUILD_TESTS

trunk/psLib/test/math/tap_psMatrix08.c

@@ -1 +1 @@
 /** @file  tap_psMatrix_08.c
 *
-*  @brief Test driver for psMatrix transpose function
-*
-*  This test driver contains the following tests:
-*     Transpose input image into output image
-*     Transpose input image into auto allocated NULL output image
-*
-*  @author  Ross Harman, MHPCC
+*  @brief psMatrixLUSolve, psMatrixGJSolve tests (for ill-conditioned matrix)
+*  @author  Eugene Magnier, IfA
 *
 *  @version $Revision: 1.2 $  $Name: not supported by cvs2svn $
 *  @date  $Date: 2007-05-02 04:20:06 $
 *
-*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+*  Copyright 2004-2005 Institute for Astronomy, University of Hawaii
 *
 */
@@ -21 +16 @@
 #include "pstap.h"
 
-# if (0)
-bool check_matrix(psImage *img)
-{
-    bool errorFlag = false;
-
-    for(psU32 i=0; i<img->numRows; i++) {
-        for(psU32 j=0; j<img->numCols; j++) {
-            if(img->type.type == PS_TYPE_F64) {
-                if(fabs(img->data.F64[i][j]-truthMatrix[i][j]) > TOLERANCE) {
-                    diag("Matrix values at element %d, %d don't agree %lf vs %lf\n", i, j,
-                         img->data.F64[i][j], truthMatrix[i][j]);
-                    errorFlag = true;
-                }
-            } else if(img->type.type == PS_TYPE_F32) {
-                if(fabs(img->data.F32[i][j]-truthMatrix[i][j]) > TOLERANCE) {
-                    diag("Matrix values at element %d, %d don't agree %f vs %lf\n", i, j,
-                         img->data.F32[i][j], truthMatrix[i][j]);
-                    errorFlag = true;
-                }
-            }
-        }
-    }
-    return(errorFlag);
-}
-# endif
+# define DEBUG 1
 
 psS32 main( psS32 argc, char* argv[] )
 {
-    plan_tests(14);
+    plan_tests(23);
+    // psTraceSetLevel("psLib.math.psMatrixGJSolve", 4);
 
     // Transpose input image into output image
@@ -55 +27 @@
         psMemId id = psMemGetId();
 
-        // we have a specific image which gave us trouble elsewhere:
-        psImage *inImage = psImageAlloc(2, 2, PS_TYPE_F64);
-        inImage->data.F64[0][0] = 0.0;
-        inImage->data.F64[0][1] = -0.000100588316855;
-        inImage->data.F64[1][0] = +0.000100588316855;
-        inImage->data.F64[1][1] = 4.4276353417e-11;
-
-        psVector *inVector = psVectorAlloc (2, PS_TYPE_F64);
-        inVector->data.F64[0] = 0.0;
-        inVector->data.F64[1] = 9.06388443347e-08;
+        psFits *fits = NULL;
+
+        // we have a specific image and vector pair which gave us trouble elsewhere:
+        // XXX this is an ill-conditioned matrix.  LU Decomposition does not inform us that it is ill-conditioned.  
+        // the result solves the equation, but what are the errors on the values?
+        fits = psFitsOpen ("data/Agj.fits", "r");
+        ok(fits, "opened test image Agj.fits");
+
+        psImage *Aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0);
+        ok(Aimage, "loaded test image Agj.fits");
+
+        psImage *aimage = psImageCopy (NULL, Aimage, Aimage->type.type);
+        ok(aimage, "copied test image Agj.fits");
+
+        psFitsClose (fits);
+
+        fits = psFitsOpen ("data/Bgj.fits", "r");
+        ok(fits, "opened test image Bgj.fits");
+
+        psImage *Bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0);
+        ok(Aimage, "loaded test image Bgj.fits");
+
+        psFitsClose (fits);
+
+        psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type);
+        ok(Bvector, "allocated B vector");
+
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = psImageGet (Bimage, 0, i);
+            psVectorSet (Bvector, i, value);
+        }
 
         bool status;
-        status = psMatrixGJSolveF32(inImage, inVector);
-        ok(status, "psMatrixGJSolve returns true exit status");
-        psFree(inImage);
-        psFree(inVector);
+        status = psMatrixLUSolve(Aimage, Bvector);
+        ok(!status, "psMatrixLUSolve correctly returns false for ill-conditioned matrix");
+
+# if (DEBUG)
+        fprintf (stderr, "LU Solution:\n");
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = psVectorGet (Bvector, i);
+            double valerr = psImageGet (Aimage, i, i);
+            fprintf (stderr, "%f +/- %f\n", value, valerr);
+        }
+
+        // calculate Ax and compare with B:
+        fprintf (stderr, "result:\n");
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = 0;
+            for (int j = 0; j < Bvector->n; j++) {
+                double tmpV = psVectorGet (Bvector, j);
+                double tmpI = psImageGet (aimage, j, i);
+                value += tmpV*tmpI;
+            }
+            double actual = psImageGet (Bimage, 0, i);
+            fprintf (stderr, "%f vs %f (delta: %f)\n", value, actual, actual - value);
+        }
+# endif
+
+        psFree(Aimage);
+        psFree(Bimage);
+        psFree(Bvector);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // Transpose input image into output image
+    {
+        psMemId id = psMemGetId();
+
+        psFits *fits = NULL;
+
+        // we have a specific ill-conditioned matrix in Agj.fits. psMatrixGJSolve detects this and reports a failure.
+        fits = psFitsOpen ("data/Agj.fits", "r");
+        ok(fits, "opened test image Agj.fits");
+
+        psImage *Aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0);
+        ok(Aimage, "loaded test image Agj.fits");
+
+        psFitsClose (fits);
+
+        fits = psFitsOpen ("data/Bgj.fits", "r");
+        ok(fits, "opened test image Bgj.fits");
+
+        psImage *Bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0);
+        ok(Bimage, "loaded test image Bgj.fits");
+
+        psFitsClose (fits);
+
+        psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type);
+        ok(Bvector, "allocated B vector");
+
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = psImageGet (Bimage, 0, i);
+            psVectorSet (Bvector, i, value);
+        }
+
+        bool status;
+        status = psMatrixGJSolve(Aimage, Bvector);
+        ok(!status, "psMatrixGJSolve correctly returns false for ill-conditioned matrix");
+
+# if (DEBUG)
+        fprintf (stderr, "GJ Solution:\n");
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = psVectorGet (Bvector, i);
+            fprintf (stderr, "%f\n", value);
+        }
+# endif
+
+        psFree(Aimage);
+        psFree(Bimage);
+        psFree(Bvector);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // Transpose input image into output image
+    {
+        psMemId id = psMemGetId();
+
+        psFits *fits = NULL;
+
+        // we have a specific ill-conditioned matrix in Agj.fits. psMatrixGJSolve detects this and reports a failure.
+        fits = psFitsOpen ("data/Agj.fits", "r");
+        ok(fits, "opened test image Agj.fits");
+
+        psImage *aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0);
+        ok(aimage, "loaded test image Agj.fits");
+
+        psImage *Aimage = psImageCopy (NULL, aimage, PS_TYPE_F64);
+        ok(Aimage, "converted test image to F64");
+
+        psFitsClose (fits);
+
+        fits = psFitsOpen ("data/Bgj.fits", "r");
+        ok(fits, "opened test image Bgj.fits");
+
+        psImage *bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0);
+        ok(bimage, "loaded test image Bgj.fits");
+
+        psImage *Bimage = psImageCopy (NULL, bimage, PS_TYPE_F64);
+        ok(Bimage, "converted test image to F64");
+
+        psFitsClose (fits);
+
+        psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type);
+        ok(Bvector, "allocated B vector");
+
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = psImageGet (Bimage, 0, i);
+            psVectorSet (Bvector, i, value);
+        }
+
+        bool status;
+        status = psMatrixGJSolve(Aimage, Bvector);
+        ok(!status, "psMatrixGJSolve correctly returns false for ill-conditioned matrix");
+
+# if (DEBUG)    
+        fprintf (stderr, "GJ Solution:\n");
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = psVectorGet (Bvector, i);
+            double valerr = psImageGet (Aimage, i, i);
+            fprintf (stderr, "%f +/- %f\n", value, valerr);
+        }
+        
+        // calculate Ax and compare with B:
+        fprintf (stderr, "result:\n");
+        for (int i = 0; i < Bimage->numRows; i++) {
+            double value = 0;
+            for (int j = 0; j < Bvector->n; j++) {
+                double tmpV = psVectorGet (Bvector, j);
+                double tmpI = psImageGet (aimage, j, i);
+                value += tmpV*tmpI;
+            }
+            double actual = psImageGet (Bimage, 0, i);
+            fprintf (stderr, "%f vs %f (delta: %f)\n", value, actual, actual - value);
+        }
+# endif
+
+        psFree(Aimage);
+        psFree(Bimage);
+        psFree(aimage);
+        psFree(bimage);
+        psFree(Bvector);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }