[27953d1] | 1 | #!/usr/bin/env python |
---|
| 2 | """ |
---|
| 3 | Class to validate a given 2D model w/ 3 axes by averaging it |
---|
| 4 | and comparing to 1D prediction. |
---|
| 5 | |
---|
| 6 | The equation used for averaging is: |
---|
| 7 | |
---|
| 8 | (integral dphi from 0 to 2pi)(integral dtheta from 0 to pi) |
---|
| 9 | p(theta, phi) I(q) sin(theta) dtheta |
---|
| 10 | |
---|
| 11 | = (1/N_phi) (1/N_theta) (pi/2) (sum over N_phi) (sum over N_theta) |
---|
| 12 | p(theta_i, phi_i) I(q) sin(theta_i) |
---|
| 13 | |
---|
| 14 | where p(theta, phi) is the probability distribution normalized to 4pi. |
---|
| 15 | In the current case, we put p(theta, phi) = 1. |
---|
| 16 | |
---|
| 17 | The normalization factor results from: |
---|
| 18 | 2pi/N_phi for the phi sum |
---|
| 19 | x pi/N_theta for the theta sum |
---|
| 20 | x 1/(4pi) because p is normalized to 4pi |
---|
| 21 | -------------- |
---|
| 22 | = (1/N_phi) (1/N_theta) (pi/2) |
---|
| 23 | |
---|
| 24 | Note: |
---|
| 25 | Averaging the 3-axes 2D scattering intensity give a slightly |
---|
| 26 | different output than the 1D function |
---|
[bda194e3] | 27 | at hight Q (Q>~0.2). This is due to the way the IGOR library |
---|
| 28 | averages(?), taking only 76 points in alpha, the angle between |
---|
[27953d1] | 29 | the axis of the ellipsoid and the q vector. |
---|
| 30 | |
---|
| 31 | Note: Core-shell and sphere models are symmetric around |
---|
| 32 | all axes and don't need to be tested in the following way. |
---|
| 33 | """ |
---|
| 34 | import sys, math |
---|
| 35 | from sans.models.EllipticalCylinderModel import EllipticalCylinderModel |
---|
| 36 | from sans.models.ParallelepipedModel import ParallelepipedModel |
---|
| 37 | from sans.models.TriaxialEllipsoidModel import TriaxialEllipsoidModel |
---|
| 38 | |
---|
| 39 | class Validate2D: |
---|
| 40 | """ |
---|
| 41 | Class to validate a given 2D model by averaging it |
---|
| 42 | and comparing to 1D prediction. |
---|
| 43 | """ |
---|
| 44 | |
---|
| 45 | def __init__(self): |
---|
| 46 | """ Initialization """ |
---|
| 47 | # Precision for the result comparison |
---|
| 48 | self.precision = 0.000001 |
---|
| 49 | # Flag for end result |
---|
| 50 | self.passed = True |
---|
| 51 | |
---|
| 52 | def __call__(self, model_class=EllipticalCylinderModel, points = 101): |
---|
| 53 | """ |
---|
| 54 | Perform test and produce output file |
---|
| 55 | @param model_class: python class of the model to test |
---|
| 56 | """ |
---|
| 57 | print "Averaging %s" % model_class.__name__ |
---|
| 58 | passed = True |
---|
| 59 | |
---|
| 60 | npts =points |
---|
[bda194e3] | 61 | #The average values are very sensitive to npts of phi so npts_alpha should be large enough. |
---|
| 62 | npts_alpha =180 #npts of phi |
---|
[27953d1] | 63 | model = model_class() |
---|
| 64 | #model.setParam('scale', 1.0) |
---|
| 65 | #model.setParam('contrast', 1.0) |
---|
| 66 | |
---|
| 67 | theta_label = 'cyl_theta' |
---|
| 68 | if not model.params.has_key(theta_label): |
---|
| 69 | theta_label = 'parallel_theta' |
---|
| 70 | if not model.params.has_key(theta_label): |
---|
| 71 | theta_label = 'axis_theta' |
---|
| 72 | |
---|
| 73 | phi_label = 'cyl_phi' |
---|
| 74 | if not model.params.has_key(phi_label): |
---|
| 75 | phi_label = 'parallel_phi' |
---|
| 76 | if not model.params.has_key(phi_label): |
---|
| 77 | phi_label = 'axis_phi' |
---|
| 78 | |
---|
| 79 | psi_label = 'cyl_psi' |
---|
| 80 | if not model.params.has_key(psi_label): |
---|
| 81 | psi_label = 'parallel_psi' |
---|
| 82 | if not model.params.has_key(psi_label): |
---|
| 83 | psi_label = 'axis_psi' |
---|
| 84 | |
---|
[8e36cdd] | 85 | output_f = open("%s_avg.txt" % model.__class__.__name__,'w') |
---|
[27953d1] | 86 | output_f.write("<q_average> <2d_average> <1d_average>\n") |
---|
| 87 | |
---|
[bda194e3] | 88 | for i_q in range(1, 23): |
---|
[27953d1] | 89 | q = 0.01*i_q |
---|
| 90 | sum = 0.0 |
---|
[bda194e3] | 91 | weight = 0.0 |
---|
[27953d1] | 92 | for i_theta in range(npts): |
---|
[4628e31] | 93 | theta = 180.0/npts*(i_theta+1) |
---|
[27953d1] | 94 | |
---|
| 95 | model.setParam(theta_label, theta) |
---|
| 96 | |
---|
[bda194e3] | 97 | for j in range(npts_alpha): |
---|
[4628e31] | 98 | model.setParam(phi_label, 180.0 * 2.0 / npts_alpha * j) |
---|
[27953d1] | 99 | for k in range(npts): |
---|
[4628e31] | 100 | model.setParam(psi_label, 180.0 * 2.0 / npts * k) |
---|
[27953d1] | 101 | if str(model.run([q, 0])).count("INF")>0: |
---|
[4628e31] | 102 | print "ERROR", q, theta, 180.0 * 2.0 / npts * k |
---|
[bda194e3] | 103 | |
---|
| 104 | # sin() is due to having not uniform bin number density wrt the q plane. |
---|
[4628e31] | 105 | sum += model.run([q, 0])*math.sin(theta*math.pi/180.0) |
---|
| 106 | weight += math.sin(theta*math.pi/180.0) |
---|
[bda194e3] | 107 | |
---|
| 108 | value = sum/weight #*math.pi/2.0 |
---|
[27953d1] | 109 | ana = model.run(q) |
---|
[bda194e3] | 110 | if q<0.3 and (value-ana)/ana>0.05: |
---|
[27953d1] | 111 | passed = False |
---|
| 112 | output_f.write("%10g %10g %10g\n" % (q, value, ana)) |
---|
| 113 | print "Q=%g: %10g %10g %10g %10g" % (q, value, ana, value-ana, value/ana) |
---|
| 114 | |
---|
| 115 | output_f.close() |
---|
| 116 | return passed |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | |
---|
| 120 | if __name__ == '__main__': |
---|
| 121 | validator = Validate2D() |
---|
[8e36cdd] | 122 | |
---|
[bda194e3] | 123 | #Note: Test one model by one model, otherwise it could crash depending on the memory. |
---|
| 124 | |
---|
[27953d1] | 125 | te_passed =validator(TriaxialEllipsoidModel, points=76) |
---|
[bda194e3] | 126 | #pp_passed = validator(ParallelepipedModel, points=76) |
---|
| 127 | #ell_passed = validator(EllipticalCylinderModel, points=76) |
---|
[27953d1] | 128 | |
---|
| 129 | print "" |
---|
| 130 | print "Model Passed" |
---|
| 131 | print "TriaxialEllipsoid %s" % te_passed |
---|
[bda194e3] | 132 | #print "ParallelepipedModel %s" % pp_passed |
---|
| 133 | #print "EllipticalCylinder %s" % ell_passed |
---|
[27953d1] | 134 | |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | |
---|
| 138 | |
---|