source: sasview/sansmodels/test/validate_2D_elliptical_cyl.py @ 162c778

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 162c778 was 18e250c, checked in by Gervaise Alina <gervyh@…>, 13 years ago

move test to sansmodels top level

  • Property mode set to 100644
File size: 7.7 KB
Line 
1#!/usr/bin/env python
2"""
3    Test for Elliptical Cylinder model.
4    Class to validate a given 2D model by averaging it
5    and comparing to 1D prediction.
6   
7"""
8import sys, math
9from sans.models.EllipticalCylinderModel import EllipticalCylinderModel
10
11
12class Validate2D:
13    """
14        Class to validate a given 2D model by averaging it
15        and comparing to 1D prediction.
16    """
17   
18    def __init__(self):
19        """ Initialization """
20        # Precision for the result comparison
21        self.precision = 0.000001
22        # Flag for end result
23        self.passed = True
24        # Verbose flag
25        self.verbose = True
26               
27    def __call__(self, npts = 101):
28        """
29            Perform test and produce output file
30            @param npts: number of points to average over
31            @return: True if the test passed, otherwise False
32        """
33        passed = True
34       
35        model = EllipticalCylinderModel()
36       
37        theta_label = 'cyl_theta'
38        if not model.params.has_key(theta_label):
39            theta_label = 'axis_theta'
40           
41        phi_label = 'cyl_phi'
42        if not model.params.has_key(phi_label):
43            phi_label = 'axis_phi'
44       
45        output_f = open("average_func.txt",'w')   
46        output_f.write("<q_average> <2d_average> <1d_average>\n")
47           
48        for i_q in range(1, 15):
49            q = 0.3/15.0*i_q
50            value = self.average_point_3D(model, q, npts)
51            ana = model.run(q)
52            if q<0.3 and (value-ana)/ana>0.05:
53                passed = False
54            output_f.write("%10g %10g %10g\n" % (q, value, ana))
55            if self.verbose:
56                print "Q=%g: %10g %10g %10g %10g" % (q, value, ana, value-ana, value/ana)
57       
58        output_f.close()
59        return passed
60   
61    def average_point_3D(self, model, q, npts):
62        """
63            Average intensity over all orientations
64            of the main cylinder axis and the
65            rotation around that axis
66           
67            @param model: model to test
68            @param q: q-value
69            @param npts: number of points to average over
70            @return: average intensity
71        """
72        sum = 0.0
73        for i_theta in range(npts):
74            theta = math.pi/npts*i_theta
75           
76            model.setParam('cyl_theta', theta)
77           
78            for j in range(npts):
79                model.setParam('cyl_phi', math.pi * 2.0 / npts * j)
80                     
81                for k in range(npts):
82                    model.setParam("cyl_psi", math.pi * 2.0 / npts * k)
83                   
84                    if str(model.run([q, 0])).count("IN")>0:                       
85                        if self.verbose:
86                            #print "ERROR", q, theta, math.pi * 2.0 / npts * j
87                            pass
88                    else:
89                        sum += math.sin(theta)*model.run([q, 0])
90         
91        value = sum/npts/npts/npts
92        return value
93
94    def checkCylinder2D(self, phi):
95        """
96            Check that the 2D scattering intensity reduces
97            to a cylinder when r_ratio = 1.0
98            @param phi: angle of the vector q on the detector
99            @return: True if the test passed, otherwise False
100        """
101        from sans.models.CylinderModel import CylinderModel
102       
103        cyl = CylinderModel()
104        cyl.setParam("cyl_theta", 90)
105        cyl.setParam("cyl_phi", 0.0)
106        cyl.setParam("radius", 20)
107        cyl.setParam("length", 400)
108        cyl.setParam("sldCyl", 2.0e-6)
109        cyl.setParam("sldSolv", 1.0e-6)
110
111        ell = EllipticalCylinderModel()
112        ell.setParam("r_ratio", 1.0)
113        ell.setParam("r_minor", 20)
114        ell.setParam("cyl_theta", 90)
115        ell.setParam("cyl_phi", 0.0)
116        ell.setParam("length", 400)
117        ell.setParam("sldCyl", 2.0e-6)
118        ell.setParam("sldSolv", 1.0e-6)
119       
120        passed = True
121        for i_q in range(1, 30):
122            q = 0.025*i_q
123            ell_val = ell.run([q, phi])
124            cyl_val = cyl.run([q, phi])
125            if self.verbose:
126                print "Q=%g    Ell=%g    Cyl=%g   R=%g" %(q, ell_val, cyl_val, ell_val/cyl_val)
127            if math.fabs(ell_val-cyl_val)/cyl_val>0.05:
128                passed= False
129               
130        return passed
131
132       
133    def checkCylinder(self, points):
134        """
135            Compare the average over all orientations
136            of the main cylinder axis for a cylinder
137            and the elliptical cylinder with r_ratio = 1
138           
139            @param points: number of points to average over
140            @return: True if the test passed, otherwise False
141        """
142        from sans.models.CylinderModel import CylinderModel
143       
144        passed = True
145       
146        npts =points
147        model = EllipticalCylinderModel()
148        model.setParam('r_ratio', 1.0)
149        model.setParam("r_minor", 20)
150        model.setParam("cyl_theta", 90)
151        model.setParam("cyl_phi", 0.0)
152        model.setParam("length", 400)
153        model.setParam("sldEll", 2.0e-6)
154        model.setParam("sldSolv", 1.0e-6)
155       
156        cyl = CylinderModel()
157        cyl.setParam("cyl_theta", 90)
158        cyl.setParam("cyl_phi", 0.0)
159        cyl.setParam("radius", 20)
160        cyl.setParam("length", 400)
161        cyl.setParam("sldCyl", 2.0e-6)
162        cyl.setParam("sldSolv", 1.0e-6)
163
164       
165        output_f = open("average_func.txt",'w')   
166        output_f.write("<q_average> <2d_average> <1d_average>\n")
167           
168        for i_q in range(1, 15):
169            q = 0.3/15.0*i_q
170            value = self.average_point_2D(model, q, npts)
171           
172            ana = cyl.run(q)
173            if q<0.3 and math.fabs(value-ana)/ana>0.05:
174                passed = False
175            output_f.write("%10g %10g %10g\n" % (q, value, ana))
176            if self.verbose:
177                print "Q=%g: %10g %10g %10g %10g" % (q, value, ana, value-ana, value/ana)
178       
179        output_f.close()
180        return passed
181       
182       
183
184    def average_point_2D(self, model, q, npts):
185        """
186            Average intensity over all orientations
187            of the main cylinder axis
188           
189            @param model: model to test
190            @param q: q-value
191            @param npts: number of points to average over
192            @return: average intensity
193        """
194        sum = 0.0
195       
196        for i_theta in range(npts):
197            theta = math.pi/npts*i_theta           
198            model.setParam('cyl_theta', theta * 180 / math.pi)
199           
200            for j in range(npts):
201                model.setParam('cyl_phi', 180 * 2.0 / npts * j)
202                     
203                if str(model.run([q, 0])).count("IN")>0:                       
204                    if self.verbose:
205                        print "ERROR", q, theta, 180 * 2.0 / npts * j
206                else:
207                    sum += math.sin(theta)*model.run([q, 0])
208       
209        value = sum/npts/npts*math.pi/2.0
210        return value
211
212   
213if __name__ == '__main__':
214    select = 0
215
216    validator = Validate2D()
217    validator.verbose = True
218   
219    print "Testing Elliptical cylinder to 5%\n"
220   
221    if select == 0 or select == 1:
222        # Check that the scat intensity reduces to a cylinder
223        # for r_ratio = 1
224        print "Comparing to Cyl I(q,theta):", validator.checkCylinder2D(1.5)
225    elif select == 0 or select == 2:
226        print "Comparing average to Cyl I(q):", validator.checkCylinder(100)
227    elif select == 0 or select == 3:
228        ellcyl_passed = validator(100)
229   
230       
231       
232       
233 
234   
Note: See TracBrowser for help on using the repository browser.