I get slightly different results when I repeat a calculation.

I've seen this problem before (it went away but has returned):

http://projects.scipy.org/pipermail/numpy-discussion/2007-January/025724.html

A unit test is attached. It contains three tests:

In test1, I construct matrices x and y and then repeatedly calculate z
= calc(x,y). The result z is the same every time. So this test passes.

In test2, I construct matrices x and y each time before calculating z
= calc(x,y). Sometimes z is slightly different. But the x's test to be
equal and so do the y's. This test fails (on Debian Lenny, Core 2 Duo,
with libatlas3gf-sse2 but not with libatlas3gf-sse).

test3 is the same as test2 but I calculate z like this: z =
calc(100*x,y) / (100 * 100). This test passes.

I get:

======================================================================
FAIL: repeatability #2
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/[snip]/test/repeat_test.py", line 73, in test_repeat_2
    self.assert_(result, msg)
AssertionError: Max difference = 2.04946e-16

----------------------------------------------------------------------

Should a unit test like this be added to numpy?
"""Unit tests for repeatability.

The problem and possible causes are discussed in the following thread:
http://projects.scipy.org/pipermail/numpy-discussion/2007-January/025724.html

If the tests do no pass, then try removing ATLAS-SSE2 from your system. But
similar problems (not tested) may occur unless you remove all versions of
ATLAS (SSE and BASE, for example).
"""

import unittest
import numpy.matlib as M
M.seterr(divide='ignore')
M.seterr(invalid='ignore')

def calc(x, y): 
    return x * (x.T * y)
    
def load():

    # x data
    x = M.zeros((3,3))
    x[0,0] =  0.00301404794991108
    x[0,1] =  0.00264742266711118
    x[0,2] = -0.00112705028731085
    x[1,0] =  0.0228605377994491
    x[1,1] =  0.00337153112741583
    x[1,2] = -0.00823674912992519
    x[2,0] =  0.00447839875836716
    x[2,1] =  0.00274880280576514
    x[2,2] = -0.00161133933606597
         
    # y data
    y = M.zeros((3,1))
    y[0,0] = 0.000885398
    y[1,0] = 0.00667193
    y[2,0] = 0.000324727
    
    return x, y    

class Test_repeat(unittest.TestCase):
    "Test repeatability"
    
    def setUp(self):
        self.nsim = 100
        self.x, self.y = load()
        
    def test_repeat_1(self):
        "repeatability #1"
        z0 = calc(self.x, self.y)
        msg = ''
        result = True
        for i in xrange(self.nsim):
            z = calc(self.x, self.y)
            if (z != z0).any():
                msg =  'Max difference = %g' % abs((z - z0)/z0).max()     
                result = False 
                break
        self.assert_(result, msg)  

    def test_repeat_2(self):
        "repeatability #2" 
        z0 = calc(self.x, self.y)
        msg = ''
        result = True
        for i in xrange(self.nsim):
            x, y = load() 
            z = calc(x, y)
            if (z != z0).any():
                msg =  'Max difference = %g' % abs((z - z0)/z0).max()       
                result = False 
                break
        self.assert_(result, msg)      

    def test_repeat_3(self):
        "repeatability #3"  
        z0 = calc(100*self.x, self.y) / (100 * 100)
        msg = ''
        result = True
        for i in xrange(self.nsim):
            x, y = load() 
            z = calc(100*x, y) / (100 * 100)
            if (z != z0).any():
                msg =  'Max difference = %g' % abs((z - z0)/z0).max()       
                result = False 
                break
        self.assert_(result, msg) 

 
def testsuite():
    unit = unittest.TestLoader().loadTestsFromTestCase
    s = []
    s.append(unit(Test_repeat))    
    return unittest.TestSuite(s)

def run():   
    suite = testsuite()
    unittest.TextTestRunner(verbosity=2).run(suite)
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Reply via email to