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

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 bd30f4a5 was 9a23253e, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

Started to add slit smearing. Added Rg and I(q=0) as outputs.

  • Property mode set to 100644
File size: 17.2 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 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_inversion(self):
179        """
180            Test an inversion for which we know the answer
181        """
182        x, y, err = load("sphere_80.txt")
183
184        # Choose the right d_max...
185        self.invertor.d_max = 160.0
186        # Set a small alpha
187        self.invertor.alpha = 1e-7
188        # Set data
189        self.invertor.x   = x
190        self.invertor.y   = y
191        self.invertor.err = err
192        # Perform inversion
193        out, cov = self.invertor.invert_optimize(10)
194        #out, cov = self.invertor.invert(10)
195        # This is a very specific case
196        # We should make sure it always passes
197        self.assertTrue(self.invertor.chi2/len(x)<200.00)
198       
199        # Check the computed P(r) with the theory
200        # for shpere of radius 80
201        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
202        y = numpy.zeros(len(x))
203        dy = numpy.zeros(len(x))
204        y_true = numpy.zeros(len(x))
205
206        sum = 0.0
207        sum_true = 0.0
208        for i in range(len(x)):
209            #y[i] = self.invertor.pr(out, x[i])
210            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
211            sum += y[i]
212            if x[i]<80.0:
213                y_true[i] = pr_theory(x[i], 80.0)
214            else:
215                y_true[i] = 0
216            sum_true += y_true[i]
217           
218        y = y/sum*self.invertor.d_max/len(x)
219        dy = dy/sum*self.invertor.d_max/len(x)
220        y_true = y_true/sum_true*self.invertor.d_max/len(x)
221       
222        chi2 = 0.0
223        for i in range(len(x)):
224            res = (y[i]-y_true[i])/dy[i]
225            chi2 += res*res
226           
227        try:
228            self.assertTrue(chi2/51.0<10.0)
229        except:
230            print "chi2 =", chi2/51.0
231            raise
232       
233    def test_lstsq(self):
234        """
235            Test an inversion for which we know the answer
236        """
237        x, y, err = load("sphere_80.txt")
238
239        # Choose the right d_max...
240        self.invertor.d_max = 160.0
241        # Set a small alpha
242        self.invertor.alpha = .005
243        # Set data
244        self.invertor.x   = x
245        self.invertor.y   = y
246        self.invertor.err = err
247        # Perform inversion
248        #out, cov = self.invertor.invert(10)
249       
250        out, cov = self.invertor.lstsq(10)
251         
252       
253        # This is a very specific case
254        # We should make sure it always passes
255        try:
256            self.assertTrue(self.invertor.chi2/len(x)<200.00)
257        except:
258            print "Chi2(I(q)) =", self.invertor.chi2/len(x)
259            raise
260       
261        # Check the computed P(r) with the theory
262        # for shpere of radius 80
263        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
264        y = numpy.zeros(len(x))
265        dy = numpy.zeros(len(x))
266        y_true = numpy.zeros(len(x))
267
268        sum = 0.0
269        sum_true = 0.0
270        for i in range(len(x)):
271            #y[i] = self.invertor.pr(out, x[i])
272            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
273            sum += y[i]
274            if x[i]<80.0:
275                y_true[i] = pr_theory(x[i], 80.0)
276            else:
277                y_true[i] = 0
278            sum_true += y_true[i]
279           
280        y = y/sum*self.invertor.d_max/len(x)
281        dy = dy/sum*self.invertor.d_max/len(x)
282        y_true = y_true/sum_true*self.invertor.d_max/len(x)
283       
284        chi2 = 0.0
285        for i in range(len(x)):
286            res = (y[i]-y_true[i])/dy[i]
287            chi2 += res*res
288           
289        try:
290            self.assertTrue(chi2/51.0<50.0)
291        except:
292            print "chi2(P(r)) =", chi2/51.0
293            raise
294       
295        # Test the number of peaks
296        self.assertEqual(self.invertor.get_peaks(out), 1)
297           
298    def test_q_zero(self):
299        """
300            Test error condition where a point has q=0
301        """
302        x, y, err = load("sphere_80.txt")
303        x[0] = 0.0
304       
305        # Choose the right d_max...
306        self.invertor.d_max = 160.0
307        # Set a small alpha
308        self.invertor.alpha = 1e-7
309        # Set data
310        def doit():
311            self.invertor.x   = x
312        self.assertRaises(ValueError, doit)
313       
314                           
315    def test_q_neg(self):
316        """
317            Test error condition where a point has q<0
318        """
319        x, y, err = load("sphere_80.txt")
320        x[0] = -0.2
321       
322        # Choose the right d_max...
323        self.invertor.d_max = 160.0
324        # Set a small alpha
325        self.invertor.alpha = 1e-7
326        # Set data
327        self.invertor.x   = x
328        self.invertor.y   = y
329        self.invertor.err = err
330        # Perform inversion
331        out, cov = self.invertor.invert(4)
332       
333        try:
334            self.assertTrue(self.invertor.chi2>0)
335        except:
336            print "Chi2 =", self.invertor.chi2
337            raise
338                           
339    def test_Iq_zero(self):
340        """
341            Test error condition where a point has q<0
342        """
343        x, y, err = load("sphere_80.txt")
344        y[0] = 0.0
345       
346        # Choose the right d_max...
347        self.invertor.d_max = 160.0
348        # Set a small alpha
349        self.invertor.alpha = 1e-7
350        # Set data
351        self.invertor.x   = x
352        self.invertor.y   = y
353        self.invertor.err = err
354        # Perform inversion
355        out, cov = self.invertor.invert(4)
356       
357        try:
358            self.assertTrue(self.invertor.chi2>0)
359        except:
360            print "Chi2 =", self.invertor.chi2
361            raise
362       
363    def no_test_time(self):
364        x, y, err = load("sphere_80.txt")
365
366        # Choose the right d_max...
367        self.invertor.d_max = 160.0
368        # Set a small alpha
369        self.invertor.alpha = 1e-7
370        # Set data
371        self.invertor.x   = x
372        self.invertor.y   = y
373        self.invertor.err = err
374   
375        # time scales like nfunc**2
376        # on a Lenovo Intel Core 2 CPU T7400 @ 2.16GHz,
377        # I get time/(nfunc)**2 = 0.022 sec
378                           
379        out, cov = self.invertor.invert(15)
380        t16 = self.invertor.elapsed
381       
382        out, cov = self.invertor.invert(30)
383        t30 = self.invertor.elapsed
384       
385        t30s = t30/30.0**2
386        self.assertTrue( (t30s-t16/16.0**2)/t30s <1.2 )
387       
388    def test_clone(self):
389        self.invertor.x = self.x_in
390        clone = self.invertor.clone()
391       
392        for i in range(len(self.x_in)):
393            self.assertEqual(self.x_in[i], clone.x[i])
394       
395    def test_save(self):
396        x, y, err = load("sphere_80.txt")
397
398        # Choose the right d_max...
399        self.invertor.d_max = 160.0
400        # Set a small alpha
401        self.invertor.alpha = .0007
402        # Set data
403        self.invertor.x   = x
404        self.invertor.y   = y
405        self.invertor.err = err
406        # Perform inversion
407       
408        out, cov = self.invertor.lstsq(10)
409       
410        # Save
411        self.invertor.to_file("test_output.txt")
412   
413    def test_load(self):
414        self.invertor.from_file("test_output.txt")
415        self.assertEqual(self.invertor.d_max, 160.0)
416        self.assertEqual(self.invertor.alpha, 0.0007)
417        self.assertEqual(self.invertor.chi2, 836.797)
418        self.assertAlmostEqual(self.invertor.pr(self.invertor.out, 10.0), 903.31577041, 4)
419       
420    def test_qmin(self):
421        self.invertor.q_min = 1.0
422        self.assertEqual(self.invertor.q_min, 1.0)
423       
424        self.invertor.q_min = None
425        self.assertEqual(self.invertor.q_min, None)
426       
427                         
428    def test_qmax(self):
429        self.invertor.q_max = 1.0
430        self.assertEqual(self.invertor.q_max, 1.0)
431       
432        self.invertor.q_max = None
433        self.assertEqual(self.invertor.q_max, None)
434
435class TestErrorConditions(unittest.TestCase):
436   
437    def setUp(self):
438        self.invertor = Invertor()
439        self.invertor.d_max = 100.0
440       
441        # Test array
442        self.ntest = 5
443        self.x_in = numpy.ones(self.ntest)
444        for i in range(self.ntest):
445            self.x_in[i] = 1.0*(i+1)
446
447    def test_negative_errs(self):
448        """
449            Test an inversion for which we know the answer
450        """
451        x, y, err = load("data_error_1.txt")
452
453        # Choose the right d_max...
454        self.invertor.d_max = 160.0
455        # Set a small alpha
456        self.invertor.alpha = .0007
457        # Set data
458        self.invertor.x   = x
459        self.invertor.y   = y
460        self.invertor.err = err
461        # Perform inversion
462       
463        out, cov = self.invertor.lstsq(10)
464         
465    def test_zero_errs(self):
466        """
467            Have zero as an error should raise an exception
468        """
469        x, y, err = load("data_error_2.txt")
470
471        # Set data
472        self.invertor.x   = x
473        self.invertor.y   = y
474        self.invertor.err = err
475        # Perform inversion
476        self.assertRaises(ValueError, self.invertor.invert, 10)
477         
478       
479    def test_invalid(self):
480        """
481            Test an inversion for which we know the answer
482        """
483        x, y, err = load("data_error_1.txt")
484
485        # Set data
486        self.invertor.x   = x
487        self.invertor.y   = y
488        err = numpy.zeros(len(x)-1)
489        self.invertor.err = err
490        # Perform inversion
491        self.assertRaises(RuntimeError, self.invertor.invert, 10)
492         
493       
494    def test_zero_q(self):
495        """
496            One of the q-values is zero.
497            An exception should be raised.
498        """
499        x, y, err = load("data_error_3.txt")
500
501        # Set data
502        self.assertRaises(ValueError, self.invertor.__setattr__, 'x', x)
503       
504    def test_zero_Iq(self):
505        """
506            One of the I(q) points has a value of zero
507            Should not complain or crash.
508        """
509        x, y, err = load("data_error_4.txt")
510
511        # Set data
512        self.invertor.x   = x
513        self.invertor.y   = y
514        self.invertor.err = err
515        # Perform inversion
516        out, cov = self.invertor.lstsq(10)
517   
518    def test_negative_q(self):
519        """
520            One q value is negative.
521            Makes not sense, but should not complain or crash.
522        """
523        x, y, err = load("data_error_5.txt")
524
525        # Set data
526        self.invertor.x   = x
527        self.invertor.y   = y
528        self.invertor.err = err
529        # Perform inversion
530        out, cov = self.invertor.lstsq(10)
531   
532    def test_negative_Iq(self):
533        """
534            One I(q) value is negative.
535            Makes not sense, but should not complain or crash.
536        """
537        x, y, err = load("data_error_6.txt")
538
539        # Set data
540        self.invertor.x   = x
541        self.invertor.y   = y
542        self.invertor.err = err
543        # Perform inversion
544        out, cov = self.invertor.lstsq(10)
545       
546    def test_nodata(self):
547        """
548             No data was loaded. Should not complain.
549        """
550        out, cov = self.invertor.lstsq(10)
551   
552
553       
554def pr_theory(r, R):
555    """
556       P(r) for a sphere
557    """
558    if r<=2*R:
559        return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
560    else:
561        return 0.0
562   
563def load(path = "sphere_60_q0_2.txt"):
564    import numpy, math, sys
565    # Read the data from the data file
566    data_x   = numpy.zeros(0)
567    data_y   = numpy.zeros(0)
568    data_err = numpy.zeros(0)
569    scale    = None
570    if not path == None:
571        input_f = open(path,'r')
572        buff    = input_f.read()
573        lines   = buff.split('\n')
574        for line in lines:
575            try:
576                toks = line.split()
577                x = float(toks[0])
578                y = float(toks[1])
579                if len(toks)>2:
580                    err = float(toks[2])
581                else:
582                    if scale==None:
583                        scale = 0.15*math.sqrt(y)
584                    err = scale*math.sqrt(y)
585                data_x = numpy.append(data_x, x)
586                data_y = numpy.append(data_y, y)
587                data_err = numpy.append(data_err, err)
588            except:
589                pass
590               
591    return data_x, data_y, data_err
592       
593if __name__ == '__main__':
594    unittest.main()
Note: See TracBrowser for help on using the repository browser.