source: sasview/pr_inversion/test/utest_invertor.py @ 383189f9

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 383189f9 was 2d06beb, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

Use simple least-square fit

  • Property mode set to 100644
File size: 10.8 KB
Line 
1"""
2    Unit tests for Invertor class
3"""
4# Disable "missing docstring" complaint
5# pylint: disable-msg=C0111
6# Disable "too many methods" complaint
7# pylint: disable-msg=R0904
8
9
10import unittest, math, numpy, pylab
11from sans.pr.invertor import Invertor
12       
13class TestBasicComponent(unittest.TestCase):
14   
15    def setUp(self):
16        self.invertor = Invertor()
17        self.invertor.d_max = 100.0
18       
19        # Test array
20        self.ntest = 5
21        self.x_in = numpy.ones(self.ntest)
22        for i in range(self.ntest):
23            self.x_in[i] = 1.0*(i+1)
24
25
26    def testset_dmax(self):
27        """
28            Set and read d_max
29        """
30        value = 15.0
31        self.invertor.d_max = value
32        self.assertEqual(self.invertor.d_max, value)
33       
34    def testset_alpha(self):
35        """
36            Set and read alpha
37        """
38        value = 15.0
39        self.invertor.alpha = value
40        self.assertEqual(self.invertor.alpha, value)
41       
42    def testset_x_1(self):
43        """
44            Setting and reading the x array the hard way
45        """
46        # Set x
47        self.invertor.x = self.x_in
48       
49        # Read it back
50        npts = self.invertor.get_nx()
51        x_out = numpy.ones(npts)
52       
53        self.invertor.get_x(x_out)
54
55        for i in range(self.ntest):
56            self.assertEqual(self.x_in[i], x_out[i])
57           
58    def testset_x_2(self):
59        """
60            Setting and reading the x array the easy way
61        """
62        # Set x
63        self.invertor.x = self.x_in
64       
65        # Read it back
66        x_out = self.invertor.x
67       
68        for i in range(self.ntest):
69            self.assertEqual(self.x_in[i], x_out[i])
70       
71    def testset_y(self):
72        """
73            Setting and reading the y array the easy way
74        """
75        # Set y
76        self.invertor.y = self.x_in
77       
78        # Read it back
79        y_out = self.invertor.y
80       
81        for i in range(self.ntest):
82            self.assertEqual(self.x_in[i], y_out[i])
83       
84    def testset_err(self):
85        """
86            Setting and reading the err array the easy way
87        """
88        # Set err
89        self.invertor.err = self.x_in
90       
91        # Read it back
92        err_out = self.invertor.err
93       
94        for i in range(self.ntest):
95            self.assertEqual(self.x_in[i], err_out[i])
96       
97    def test_iq(self):
98        """
99            Test iq calculation
100        """
101        q = 0.11
102        v1 = 8.0*math.pi**2/q * self.invertor.d_max *math.sin(q*self.invertor.d_max)
103        v1 /= ( math.pi**2 - (q*self.invertor.d_max)**2.0 )
104       
105        pars = numpy.ones(1)
106        self.assertAlmostEqual(self.invertor.iq(pars, q), v1, 2)
107       
108    def test_pr(self):
109        """
110            Test pr calculation
111        """
112        r = 10.0
113        v1 = 2.0*r*math.sin(math.pi*r/self.invertor.d_max)
114        pars = numpy.ones(1)
115        self.assertAlmostEqual(self.invertor.pr(pars, r), v1, 2)
116       
117    def test_getsetters(self):
118        self.invertor.new_data = 1.0
119        self.assertEqual(self.invertor.new_data, 1.0)
120       
121        self.assertEqual(self.invertor.test_no_data, None)
122       
123    def test_inversion(self):
124        """
125            Test an inversion for which we know the answer
126        """
127        x, y, err = load("sphere_80.txt")
128
129        # Choose the right d_max...
130        self.invertor.d_max = 160.0
131        # Set a small alpha
132        self.invertor.alpha = 1e-7
133        # Set data
134        self.invertor.x   = x
135        self.invertor.y   = y
136        self.invertor.err = err
137        # Perform inversion
138        out, cov = self.invertor.invert(10)
139        # This is a very specific case
140        # We should make sure it always passes
141        self.assertTrue(self.invertor.chi2/len(x)<200.00)
142       
143        # Check the computed P(r) with the theory
144        # for shpere of radius 80
145        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
146        y = numpy.zeros(len(x))
147        dy = numpy.zeros(len(x))
148        y_true = numpy.zeros(len(x))
149
150        sum = 0.0
151        sum_true = 0.0
152        for i in range(len(x)):
153            #y[i] = self.invertor.pr(out, x[i])
154            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
155            sum += y[i]
156            if x[i]<80.0:
157                y_true[i] = pr_theory(x[i], 80.0)
158            else:
159                y_true[i] = 0
160            sum_true += y_true[i]
161           
162        y = y/sum*self.invertor.d_max/len(x)
163        dy = dy/sum*self.invertor.d_max/len(x)
164        y_true = y_true/sum_true*self.invertor.d_max/len(x)
165       
166        chi2 = 0.0
167        for i in range(len(x)):
168            res = (y[i]-y_true[i])/dy[i]
169            chi2 += res*res
170           
171        try:
172            self.assertTrue(chi2/51.0<10.0)
173        except:
174            print "chi2 =", chi2/51.0
175            raise
176       
177    def test_lstsq(self):
178        """
179            Test an inversion for which we know the answer
180        """
181        x, y, err = load("sphere_80.txt")
182
183        # Choose the right d_max...
184        self.invertor.d_max = 160.0
185        # Set a small alpha
186        self.invertor.alpha = .0007
187        # Set data
188        self.invertor.x   = x
189        self.invertor.y   = y
190        self.invertor.err = err
191        # Perform inversion
192        #out, cov = self.invertor.invert(10)
193       
194        out, cov = self.invertor.lstsq(10)
195         
196       
197        # This is a very specific case
198        # We should make sure it always passes
199        try:
200            self.assertTrue(self.invertor.chi2/len(x)<200.00)
201        except:
202            print "Chi2(I(q)) =", self.invertor.chi2/len(x)
203            raise
204       
205        # Check the computed P(r) with the theory
206        # for shpere of radius 80
207        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
208        y = numpy.zeros(len(x))
209        dy = numpy.zeros(len(x))
210        y_true = numpy.zeros(len(x))
211
212        sum = 0.0
213        sum_true = 0.0
214        for i in range(len(x)):
215            #y[i] = self.invertor.pr(out, x[i])
216            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
217            sum += y[i]
218            if x[i]<80.0:
219                y_true[i] = pr_theory(x[i], 80.0)
220            else:
221                y_true[i] = 0
222            sum_true += y_true[i]
223           
224        y = y/sum*self.invertor.d_max/len(x)
225        dy = dy/sum*self.invertor.d_max/len(x)
226        y_true = y_true/sum_true*self.invertor.d_max/len(x)
227       
228        chi2 = 0.0
229        for i in range(len(x)):
230            res = (y[i]-y_true[i])/dy[i]
231            chi2 += res*res
232           
233        try:
234            self.assertTrue(chi2/51.0<50.0)
235        except:
236            print "chi2(P(r)) =", chi2/51.0
237            raise
238           
239    def test_q_zero(self):
240        """
241            Test error condition where a point has q=0
242        """
243        x, y, err = load("sphere_80.txt")
244        x[0] = 0.0
245       
246        # Choose the right d_max...
247        self.invertor.d_max = 160.0
248        # Set a small alpha
249        self.invertor.alpha = 1e-7
250        # Set data
251        def doit():
252            self.invertor.x   = x
253        self.assertRaises(ValueError, doit)
254       
255                           
256    def test_q_neg(self):
257        """
258            Test error condition where a point has q<0
259        """
260        x, y, err = load("sphere_80.txt")
261        x[0] = -0.2
262       
263        # Choose the right d_max...
264        self.invertor.d_max = 160.0
265        # Set a small alpha
266        self.invertor.alpha = 1e-7
267        # Set data
268        self.invertor.x   = x
269        self.invertor.y   = y
270        self.invertor.err = err
271        # Perform inversion
272        out, cov = self.invertor.invert(4)
273       
274        try:
275            self.assertTrue(self.invertor.chi2>0)
276        except:
277            print "Chi2 =", self.invertor.chi2
278            raise
279                           
280    def test_Iq_zero(self):
281        """
282            Test error condition where a point has q<0
283        """
284        x, y, err = load("sphere_80.txt")
285        y[0] = 0.0
286       
287        # Choose the right d_max...
288        self.invertor.d_max = 160.0
289        # Set a small alpha
290        self.invertor.alpha = 1e-7
291        # Set data
292        self.invertor.x   = x
293        self.invertor.y   = y
294        self.invertor.err = err
295        # Perform inversion
296        out, cov = self.invertor.invert(4)
297       
298        try:
299            self.assertTrue(self.invertor.chi2>0)
300        except:
301            print "Chi2 =", self.invertor.chi2
302            raise
303       
304    def no_test_time(self):
305        x, y, err = load("sphere_80.txt")
306
307        # Choose the right d_max...
308        self.invertor.d_max = 160.0
309        # Set a small alpha
310        self.invertor.alpha = 1e-7
311        # Set data
312        self.invertor.x   = x
313        self.invertor.y   = y
314        self.invertor.err = err
315   
316        # time scales like nfunc**2
317        # on a Lenovo Intel Core 2 CPU T7400 @ 2.16GHz,
318        # I get time/(nfunc)**2 = 0.022 sec
319                           
320        out, cov = self.invertor.invert(15)
321        t16 = self.invertor.elapsed
322       
323        out, cov = self.invertor.invert(30)
324        t30 = self.invertor.elapsed
325       
326        t30s = t30/30.0**2
327        self.assertTrue( (t30s-t16/16.0**2)/t30s <1.2 )
328       
329    def test_clone(self):
330        self.invertor.x = self.x_in
331        clone = self.invertor.clone()
332       
333        for i in range(len(self.x_in)):
334            self.assertEqual(self.x_in[i], clone.x[i])
335       
336def pr_theory(r, R):
337    """
338       P(r) for a sphere
339    """
340    if r<=2*R:
341        return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
342    else:
343        return 0.0
344   
345def load(path = "sphere_60_q0_2.txt"):
346    import numpy, math, sys
347    # Read the data from the data file
348    data_x   = numpy.zeros(0)
349    data_y   = numpy.zeros(0)
350    data_err = numpy.zeros(0)
351    if not path == None:
352        input_f = open(path,'r')
353        buff    = input_f.read()
354        lines   = buff.split('\n')
355        for line in lines:
356            try:
357                toks = line.split()
358                x = float(toks[0])
359                y = float(toks[1])
360                data_x = numpy.append(data_x, x)
361                data_y = numpy.append(data_y, y)
362                # Set the error of the first point to 5%
363                # to make theory look like data
364                scale = 0.1/math.sqrt(data_x[0])
365                data_err = numpy.append(data_err, scale*math.sqrt(y))
366            except:
367                pass
368               
369    return data_x, data_y, data_err
370       
371if __name__ == '__main__':
372    unittest.main()
Note: See TracBrowser for help on using the repository browser.