source: sasview/test/pr_inversion/test/utest_invertor.py @ 3a473ef

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 3a473ef was 3a473ef, checked in by Piotr Rozyczko <rozyczko@…>, 14 months ago

Manual update of test/ directory with changes on master - SASVIEW-996

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