[ae3ce4e] | 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 | """ |
---|
| 8 | import sys, math |
---|
| 9 | from sans.models.EllipticalCylinderModel import EllipticalCylinderModel |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | class 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", 1.57) |
---|
| 105 | cyl.setParam("cyl_phi", 0.0) |
---|
| 106 | cyl.setParam("radius", 20) |
---|
| 107 | cyl.setParam("length", 400) |
---|
| 108 | cyl.setParam("contrast", 1.0e-6) |
---|
| 109 | |
---|
| 110 | ell = EllipticalCylinderModel() |
---|
| 111 | ell.setParam("r_ratio", 1.0) |
---|
| 112 | ell.setParam("r_minor", 20) |
---|
| 113 | ell.setParam("cyl_theta", 1.57) |
---|
| 114 | ell.setParam("cyl_phi", 0.0) |
---|
| 115 | ell.setParam("length", 400) |
---|
| 116 | ell.setParam("contrast", 1.0e-6) |
---|
| 117 | |
---|
| 118 | passed = True |
---|
| 119 | for i_q in range(1, 30): |
---|
| 120 | q = 0.025*i_q |
---|
| 121 | ell_val = ell.run([q, phi]) |
---|
| 122 | cyl_val = cyl.run([q, phi]) |
---|
| 123 | if self.verbose: |
---|
| 124 | print "Q=%g Ell=%g Cyl=%g R=%g" %(q, ell_val, cyl_val, ell_val/cyl_val) |
---|
| 125 | if math.fabs(ell_val-cyl_val)/cyl_val>0.05: |
---|
| 126 | passed= False |
---|
| 127 | |
---|
| 128 | return passed |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | def checkCylinder(self, points): |
---|
| 132 | """ |
---|
| 133 | Compare the average over all orientations |
---|
| 134 | of the main cylinder axis for a cylinder |
---|
| 135 | and the elliptical cylinder with r_ratio = 1 |
---|
| 136 | |
---|
| 137 | @param points: number of points to average over |
---|
| 138 | @return: True if the test passed, otherwise False |
---|
| 139 | """ |
---|
| 140 | from sans.models.CylinderModel import CylinderModel |
---|
| 141 | |
---|
| 142 | passed = True |
---|
| 143 | |
---|
| 144 | npts =points |
---|
| 145 | model = EllipticalCylinderModel() |
---|
| 146 | model.setParam('r_ratio', 1.0) |
---|
| 147 | model.setParam("r_minor", 20) |
---|
| 148 | model.setParam("cyl_theta", 1.57) |
---|
| 149 | model.setParam("cyl_phi", 0.0) |
---|
| 150 | model.setParam("length", 400) |
---|
| 151 | model.setParam("contrast", 1.0e-6) |
---|
| 152 | |
---|
| 153 | cyl = CylinderModel() |
---|
| 154 | cyl.setParam("cyl_theta", 1.57) |
---|
| 155 | cyl.setParam("cyl_phi", 0.0) |
---|
| 156 | cyl.setParam("radius", 20) |
---|
| 157 | cyl.setParam("length", 400) |
---|
| 158 | cyl.setParam("contrast", 1.0e-6) |
---|
| 159 | |
---|
| 160 | output_f = open("average_func.txt",'w') |
---|
| 161 | output_f.write("<q_average> <2d_average> <1d_average>\n") |
---|
| 162 | |
---|
| 163 | for i_q in range(1, 15): |
---|
| 164 | q = 0.3/15.0*i_q |
---|
| 165 | value = self.average_point_2D(model, q, npts) |
---|
| 166 | |
---|
| 167 | ana = cyl.run(q) |
---|
| 168 | if q<0.3 and math.fabs(value-ana)/ana>0.05: |
---|
| 169 | passed = False |
---|
| 170 | output_f.write("%10g %10g %10g\n" % (q, value, ana)) |
---|
| 171 | if self.verbose: |
---|
| 172 | print "Q=%g: %10g %10g %10g %10g" % (q, value, ana, value-ana, value/ana) |
---|
| 173 | |
---|
| 174 | output_f.close() |
---|
| 175 | return passed |
---|
| 176 | |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | def average_point_2D(self, model, q, npts): |
---|
| 180 | """ |
---|
| 181 | Average intensity over all orientations |
---|
| 182 | of the main cylinder axis |
---|
| 183 | |
---|
| 184 | @param model: model to test |
---|
| 185 | @param q: q-value |
---|
| 186 | @param npts: number of points to average over |
---|
| 187 | @return: average intensity |
---|
| 188 | """ |
---|
| 189 | sum = 0.0 |
---|
| 190 | |
---|
| 191 | for i_theta in range(npts): |
---|
| 192 | theta = math.pi/npts*i_theta |
---|
| 193 | model.setParam('cyl_theta', theta) |
---|
| 194 | |
---|
| 195 | for j in range(npts): |
---|
| 196 | model.setParam('cyl_phi', math.pi * 2.0 / npts * j) |
---|
| 197 | |
---|
| 198 | if str(model.run([q, 0])).count("IN")>0: |
---|
| 199 | if self.verbose: |
---|
| 200 | print "ERROR", q, theta, math.pi * 2.0 / npts * j |
---|
| 201 | else: |
---|
| 202 | sum += math.sin(theta)*model.run([q, 0]) |
---|
| 203 | |
---|
| 204 | value = sum/npts/npts*math.pi/2.0 |
---|
| 205 | return value |
---|
| 206 | |
---|
| 207 | |
---|
| 208 | if __name__ == '__main__': |
---|
| 209 | select = 0 |
---|
| 210 | |
---|
| 211 | validator = Validate2D() |
---|
| 212 | validator.verbose = True |
---|
| 213 | |
---|
| 214 | print "Testing Elliptical cylinder to 5%\n" |
---|
| 215 | |
---|
| 216 | if select == 0 or select == 1: |
---|
| 217 | # Check that the scat intensity reduces to a cylinder |
---|
| 218 | # for r_ratio = 1 |
---|
| 219 | print "Comparing to Cyl I(q,theta):", validator.checkCylinder2D(1.5) |
---|
| 220 | elif select == 0 or select == 2: |
---|
| 221 | print "Comparing average to Cyl I(q):", validator.checkCylinder(100) |
---|
| 222 | elif select == 0 or select == 3: |
---|
| 223 | ellcyl_passed = validator(100) |
---|
| 224 | |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | |
---|
| 228 | |
---|
| 229 | |
---|