source: sasview/pr_inversion/test/utest_invertor.py @ ea59943

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 ea59943 was ea59943, checked in by Mathieu Doucet <doucetm@…>, 12 years ago

Re #5 improving tests

  • Property mode set to 100644
File size: 17.5 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
11from sans.pr.invertor import Invertor
12       
13class TestFiguresOfMerit(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        x, y, err = load("sphere_80.txt")
26
27        # Choose the right d_max...
28        self.invertor.d_max = 160.0
29        # Set a small alpha
30        self.invertor.alpha = .0007
31        # Set data
32        self.invertor.x   = x
33        self.invertor.y   = y
34        self.invertor.err = err
35        # Perform inversion
36        #out, cov = self.invertor.invert(10)
37       
38        self.out, self.cov = self.invertor.lstsq(10)
39
40    def test_positive(self):
41        """
42            Test whether P(r) is positive
43        """
44        self.assertEqual(self.invertor.get_positive(self.out), 1)
45       
46    def test_positive_err(self):
47        """
48            Test whether P(r) is at least 1 sigma greater than zero
49            for all r-values
50        """
51        self.assertTrue(self.invertor.get_pos_err(self.out, self.cov)>0.9)
52       
53class TestBasicComponent(unittest.TestCase):
54   
55    def setUp(self):
56        self.invertor = Invertor()
57        self.invertor.d_max = 100.0
58       
59        # Test array
60        self.ntest = 5
61        self.x_in = numpy.ones(self.ntest)
62        for i in range(self.ntest):
63            self.x_in[i] = 1.0*(i+1)
64
65    def test_has_bck_flag(self):
66        """
67            Tests the has_bck flag operations
68        """
69        self.assertEqual(self.invertor.has_bck, False)
70        self.invertor.has_bck=True
71        self.assertEqual(self.invertor.has_bck, True)
72        def doit_float():
73            self.invertor.has_bck  = 2.0
74        def doit_str():
75            self.invertor.has_bck  = 'a'
76
77        self.assertRaises(ValueError, doit_float)
78        self.assertRaises(ValueError, doit_str)
79       
80
81    def testset_dmax(self):
82        """
83            Set and read d_max
84        """
85        value = 15.0
86        self.invertor.d_max = value
87        self.assertEqual(self.invertor.d_max, value)
88       
89    def testset_alpha(self):
90        """
91            Set and read alpha
92        """
93        value = 15.0
94        self.invertor.alpha = value
95        self.assertEqual(self.invertor.alpha, value)
96       
97    def testset_x_1(self):
98        """
99            Setting and reading the x array the hard way
100        """
101        # Set x
102        self.invertor.x = self.x_in
103       
104        # Read it back
105        npts = self.invertor.get_nx()
106        x_out = numpy.ones(npts)
107       
108        self.invertor.get_x(x_out)
109
110        for i in range(self.ntest):
111            self.assertEqual(self.x_in[i], x_out[i])
112           
113    def testset_x_2(self):
114        """
115            Setting and reading the x array the easy way
116        """
117        # Set x
118        self.invertor.x = self.x_in
119       
120        # Read it back
121        x_out = self.invertor.x
122       
123        for i in range(self.ntest):
124            self.assertEqual(self.x_in[i], x_out[i])
125       
126    def testset_y(self):
127        """
128            Setting and reading the y array the easy way
129        """
130        # Set y
131        self.invertor.y = self.x_in
132       
133        # Read it back
134        y_out = self.invertor.y
135       
136        for i in range(self.ntest):
137            self.assertEqual(self.x_in[i], y_out[i])
138       
139    def testset_err(self):
140        """
141            Setting and reading the err array the easy way
142        """
143        # Set err
144        self.invertor.err = self.x_in
145       
146        # Read it back
147        err_out = self.invertor.err
148       
149        for i in range(self.ntest):
150            self.assertEqual(self.x_in[i], err_out[i])
151       
152    def test_iq(self):
153        """
154            Test iq calculation
155        """
156        q = 0.11
157        v1 = 8.0*math.pi**2/q * self.invertor.d_max *math.sin(q*self.invertor.d_max)
158        v1 /= ( math.pi**2 - (q*self.invertor.d_max)**2.0 )
159       
160        pars = numpy.ones(1)
161        self.assertAlmostEqual(self.invertor.iq(pars, q), v1, 2)
162       
163    def test_pr(self):
164        """
165            Test pr calculation
166        """
167        r = 10.0
168        v1 = 2.0*r*math.sin(math.pi*r/self.invertor.d_max)
169        pars = numpy.ones(1)
170        self.assertAlmostEqual(self.invertor.pr(pars, r), v1, 2)
171       
172    def test_getsetters(self):
173        self.invertor.new_data = 1.0
174        self.assertEqual(self.invertor.new_data, 1.0)
175       
176        self.assertEqual(self.invertor.test_no_data, None)
177       
178    def test_slitsettings(self):
179        self.invertor.slit_width = 1.0
180        self.assertEqual(self.invertor.slit_width, 1.0)
181        self.invertor.slit_height = 2.0
182        self.assertEqual(self.invertor.slit_height, 2.0)
183       
184       
185    def test_inversion(self):
186        """
187            Test an inversion for which we know the answer
188        """
189        x, y, err = load("sphere_80.txt")
190
191        # Choose the right d_max...
192        self.invertor.d_max = 160.0
193        # Set a small alpha
194        self.invertor.alpha = 1e-7
195        # Set data
196        self.invertor.x   = x
197        self.invertor.y   = y
198        self.invertor.err = err
199        # Perform inversion
200        out, cov = self.invertor.invert_optimize(10)
201        #out, cov = self.invertor.invert(10)
202        # This is a very specific case
203        # We should make sure it always passes
204        self.assertTrue(self.invertor.chi2/len(x)<200.00)
205       
206        # Check the computed P(r) with the theory
207        # for shpere of radius 80
208        x = numpy.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
209        y = numpy.zeros(len(x))
210        dy = numpy.zeros(len(x))
211        y_true = numpy.zeros(len(x))
212
213        sum = 0.0
214        sum_true = 0.0
215        for i in range(len(x)):
216            #y[i] = self.invertor.pr(out, x[i])
217            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
218            sum += y[i]
219            if x[i]<80.0:
220                y_true[i] = pr_theory(x[i], 80.0)
221            else:
222                y_true[i] = 0
223            sum_true += y_true[i]
224           
225        y = y/sum*self.invertor.d_max/len(x)
226        dy = dy/sum*self.invertor.d_max/len(x)
227        y_true = y_true/sum_true*self.invertor.d_max/len(x)
228       
229        chi2 = 0.0
230        for i in range(len(x)):
231            res = (y[i]-y_true[i])/dy[i]
232            chi2 += res*res
233           
234        try:
235            self.assertTrue(chi2/51.0<10.0)
236        except:
237            print "chi2 =", chi2/51.0
238            raise
239       
240    def test_lstsq(self):
241        """
242            Test an inversion for which we know the answer
243        """
244        x, y, err = load("sphere_80.txt")
245
246        # Choose the right d_max...
247        self.invertor.d_max = 160.0
248        # Set a small alpha
249        self.invertor.alpha = .005
250        # Set data
251        self.invertor.x   = x
252        self.invertor.y   = y
253        self.invertor.err = err
254        # Perform inversion
255        #out, cov = self.invertor.invert(10)
256       
257        out, cov = self.invertor.lstsq(10)
258         
259       
260        # This is a very specific case
261        # We should make sure it always passes
262        try:
263            self.assertTrue(self.invertor.chi2/len(x)<200.00)
264        except:
265            print "Chi2(I(q)) =", self.invertor.chi2/len(x)
266            raise
267       
268        # Check the computed P(r) with the theory
269        # for shpere of radius 80
270        x = numpy.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
271        y = numpy.zeros(len(x))
272        dy = numpy.zeros(len(x))
273        y_true = numpy.zeros(len(x))
274
275        sum = 0.0
276        sum_true = 0.0
277        for i in range(len(x)):
278            #y[i] = self.invertor.pr(out, x[i])
279            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
280            sum += y[i]
281            if x[i]<80.0:
282                y_true[i] = pr_theory(x[i], 80.0)
283            else:
284                y_true[i] = 0
285            sum_true += y_true[i]
286           
287        y = y/sum*self.invertor.d_max/len(x)
288        dy = dy/sum*self.invertor.d_max/len(x)
289        y_true = y_true/sum_true*self.invertor.d_max/len(x)
290       
291        chi2 = 0.0
292        for i in range(len(x)):
293            res = (y[i]-y_true[i])/dy[i]
294            chi2 += res*res
295           
296        try:
297            self.assertTrue(chi2/51.0<50.0)
298        except:
299            print "chi2(P(r)) =", chi2/51.0
300            raise
301       
302        # Test the number of peaks
303        self.assertEqual(self.invertor.get_peaks(out), 1)
304           
305    def test_q_zero(self):
306        """
307            Test error condition where a point has q=0
308        """
309        x, y, err = load("sphere_80.txt")
310        x[0] = 0.0
311       
312        # Choose the right d_max...
313        self.invertor.d_max = 160.0
314        # Set a small alpha
315        self.invertor.alpha = 1e-7
316        # Set data
317        def doit():
318            self.invertor.x   = x
319        self.assertRaises(ValueError, doit)
320       
321                           
322    def test_q_neg(self):
323        """
324            Test error condition where a point has q<0
325        """
326        x, y, err = load("sphere_80.txt")
327        x[0] = -0.2
328       
329        # Choose the right d_max...
330        self.invertor.d_max = 160.0
331        # Set a small alpha
332        self.invertor.alpha = 1e-7
333        # Set data
334        self.invertor.x   = x
335        self.invertor.y   = y
336        self.invertor.err = err
337        # Perform inversion
338        out, cov = self.invertor.invert(4)
339       
340        try:
341            self.assertTrue(self.invertor.chi2>0)
342        except:
343            print "Chi2 =", self.invertor.chi2
344            raise
345                           
346    def test_Iq_zero(self):
347        """
348            Test error condition where a point has q<0
349        """
350        x, y, err = load("sphere_80.txt")
351        y[0] = 0.0
352       
353        # Choose the right d_max...
354        self.invertor.d_max = 160.0
355        # Set a small alpha
356        self.invertor.alpha = 1e-7
357        # Set data
358        self.invertor.x   = x
359        self.invertor.y   = y
360        self.invertor.err = err
361        # Perform inversion
362        out, cov = self.invertor.invert(4)
363       
364        try:
365            self.assertTrue(self.invertor.chi2>0)
366        except:
367            print "Chi2 =", self.invertor.chi2
368            raise
369       
370    def no_test_time(self):
371        x, y, err = load("sphere_80.txt")
372
373        # Choose the right d_max...
374        self.invertor.d_max = 160.0
375        # Set a small alpha
376        self.invertor.alpha = 1e-7
377        # Set data
378        self.invertor.x   = x
379        self.invertor.y   = y
380        self.invertor.err = err
381   
382        # time scales like nfunc**2
383        # on a Lenovo Intel Core 2 CPU T7400 @ 2.16GHz,
384        # I get time/(nfunc)**2 = 0.022 sec
385                           
386        out, cov = self.invertor.invert(15)
387        t16 = self.invertor.elapsed
388       
389        out, cov = self.invertor.invert(30)
390        t30 = self.invertor.elapsed
391       
392        t30s = t30/30.0**2
393        self.assertTrue( (t30s-t16/16.0**2)/t30s <1.2 )
394       
395    def test_clone(self):
396        self.invertor.x = self.x_in
397        clone = self.invertor.clone()
398       
399        for i in range(len(self.x_in)):
400            self.assertEqual(self.x_in[i], clone.x[i])
401       
402    def test_save(self):
403        x, y, err = load("sphere_80.txt")
404
405        # Choose the right d_max...
406        self.invertor.d_max = 160.0
407        # Set a small alpha
408        self.invertor.alpha = .0007
409        # Set data
410        self.invertor.x   = x
411        self.invertor.y   = y
412        self.invertor.err = err
413        # Perform inversion
414       
415        out, cov = self.invertor.lstsq(10)
416       
417        # Save
418        self.invertor.to_file("test_output.txt")
419   
420    def test_load(self):
421        self.invertor.from_file("test_output.txt")
422        self.assertEqual(self.invertor.d_max, 160.0)
423        self.assertEqual(self.invertor.alpha, 0.0007)
424        self.assertEqual(self.invertor.chi2, 836.797)
425        self.assertAlmostEqual(self.invertor.pr(self.invertor.out, 10.0), 903.31577041, 4)
426       
427    def test_qmin(self):
428        self.invertor.q_min = 1.0
429        self.assertEqual(self.invertor.q_min, 1.0)
430       
431        self.invertor.q_min = None
432        self.assertEqual(self.invertor.q_min, None)
433       
434                         
435    def test_qmax(self):
436        self.invertor.q_max = 1.0
437        self.assertEqual(self.invertor.q_max, 1.0)
438       
439        self.invertor.q_max = None
440        self.assertEqual(self.invertor.q_max, None)
441
442class TestErrorConditions(unittest.TestCase):
443   
444    def setUp(self):
445        self.invertor = Invertor()
446        self.invertor.d_max = 100.0
447       
448        # Test array
449        self.ntest = 5
450        self.x_in = numpy.ones(self.ntest)
451        for i in range(self.ntest):
452            self.x_in[i] = 1.0*(i+1)
453
454    def test_negative_errs(self):
455        """
456            Test an inversion for which we know the answer
457        """
458        x, y, err = load("data_error_1.txt")
459
460        # Choose the right d_max...
461        self.invertor.d_max = 160.0
462        # Set a small alpha
463        self.invertor.alpha = .0007
464        # Set data
465        self.invertor.x   = x
466        self.invertor.y   = y
467        self.invertor.err = err
468        # Perform inversion
469       
470        out, cov = self.invertor.lstsq(10)
471         
472    def test_zero_errs(self):
473        """
474            Have zero as an error should raise an exception
475        """
476        x, y, err = load("data_error_2.txt")
477
478        # Set data
479        self.invertor.x   = x
480        self.invertor.y   = y
481        self.invertor.err = err
482        # Perform inversion
483        self.assertRaises(RuntimeError, self.invertor.invert, 10)
484       
485    def test_invalid(self):
486        """
487            Test an inversion for which we know the answer
488        """
489        x, y, err = load("data_error_1.txt")
490
491        # Set data
492        self.invertor.x   = x
493        self.invertor.y   = y
494        err = numpy.zeros(len(x)-1)
495        self.invertor.err = err
496        # Perform inversion
497        self.assertRaises(RuntimeError, self.invertor.invert, 10)
498         
499       
500    def test_zero_q(self):
501        """
502            One of the q-values is zero.
503            An exception should be raised.
504        """
505        x, y, err = load("data_error_3.txt")
506
507        # Set data
508        self.assertRaises(ValueError, self.invertor.__setattr__, 'x', x)
509       
510    def test_zero_Iq(self):
511        """
512            One of the I(q) points has a value of zero
513            Should not complain or crash.
514        """
515        x, y, err = load("data_error_4.txt")
516
517        # Set data
518        self.invertor.x   = x
519        self.invertor.y   = y
520        self.invertor.err = err
521        # Perform inversion
522        out, cov = self.invertor.lstsq(10)
523   
524    def test_negative_q(self):
525        """
526            One q value is negative.
527            Makes not sense, but should not complain or crash.
528        """
529        x, y, err = load("data_error_5.txt")
530
531        # Set data
532        self.invertor.x   = x
533        self.invertor.y   = y
534        self.invertor.err = err
535        # Perform inversion
536        out, cov = self.invertor.lstsq(10)
537   
538    def test_negative_Iq(self):
539        """
540            One I(q) value is negative.
541            Makes not sense, but should not complain or crash.
542        """
543        x, y, err = load("data_error_6.txt")
544
545        # Set data
546        self.invertor.x   = x
547        self.invertor.y   = y
548        self.invertor.err = err
549        # Perform inversion
550        out, cov = self.invertor.lstsq(10)
551       
552    def test_nodata(self):
553        """
554             No data was loaded. Should not complain.
555        """
556        out, cov = self.invertor.lstsq(10)
557   
558
559       
560def pr_theory(r, R):
561    """
562       P(r) for a sphere
563    """
564    if r<=2*R:
565        return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
566    else:
567        return 0.0
568   
569def load(path = "sphere_60_q0_2.txt"):
570    import numpy, math, sys
571    # Read the data from the data file
572    data_x   = numpy.zeros(0)
573    data_y   = numpy.zeros(0)
574    data_err = numpy.zeros(0)
575    scale    = None
576    if not path == None:
577        input_f = open(path,'r')
578        buff    = input_f.read()
579        lines   = buff.split('\n')
580        for line in lines:
581            try:
582                toks = line.split()
583                x = float(toks[0])
584                y = float(toks[1])
585                if len(toks)>2:
586                    err = float(toks[2])
587                else:
588                    if scale==None:
589                        scale = 0.15*math.sqrt(y)
590                    err = scale*math.sqrt(y)
591                data_x = numpy.append(data_x, x)
592                data_y = numpy.append(data_y, y)
593                data_err = numpy.append(data_err, err)
594            except:
595                pass
596               
597    return data_x, data_y, data_err
598       
599if __name__ == '__main__':
600    unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))
Note: See TracBrowser for help on using the repository browser.