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

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.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since cb62bd5 was cb62bd5, checked in by lewis, 2 years ago

Use manually inputted background level in Pr calculation

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