Changeset a2c1196 in sasview for sansrealspace/src

Timestamp:

Nov 2, 2007 11:02:35 AM (17 years ago)

Author:

Mathieu Doucet <doucetm@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

7e845ea

Parents:

8c050c1

Message:

Added 2D (oriented systems) error estimation, tests and validation.

Location:

sansrealspace/src/realspace

Files:

• VolumeCanvas.py (modified) (4 diffs)
• test/coreshell_d=0.01_Iq.txt (added)
• test/coreshell_d=0.01_Iq2D.txt (added)
• test/coreshell_d=0.1_Iq2D.txt (added)
• test/coreshell_d=10_Iq2D.txt (added)
• test/coreshell_d=20_Iq2D.txt (added)
• test/coreshell_d=5_Iq2D.txt (added)
• test/cylinder_d=0.5_Iq2D.txt (added)
• test/cylinder_d=1_Iq2D.txt (added)
• test/early_test.py (modified) (1 diff)
• test/ellipsoid_d=0.5_Iq2D.txt (added)
• test/sim_validation.py (modified) (5 diffs)
• test/sphere_d=2_Iq2D.txt (added)
• test/utest_oriented.py (modified) (7 diffs)

sansrealspace/src/realspace/VolumeCanvas.py

@@ -551 +551 @@
             raise ValueError, "runXY(q): bad type for q"
     
-        
-    def getIq2D(self, qx, qy):
+    def _create_modelObject(self):
         """
             Returns simulate I(q) for given q_x and q_y values.
+            Also returns model object
             @param qx: q_x [A-1]
             @param qy: q_y [A-1]
-            @return: I(q) [cm-1]
+            @return: I(q) [cm-1], model object
         """
         # To find a complete example of the correct call order:
@@ -581 +581 @@
         #pointsmodelpy.get_lorespoints(self.lores_model, self.points)
         self.npts = pointsmodelpy.get_complexpoints(self.complex_model, self.points)
+        
+        
+    def getIq2D(self, qx, qy):
+        """
+            Returns simulate I(q) for given q_x and q_y values.
+            @param qx: q_x [A-1]
+            @param qy: q_y [A-1]
+            @return: I(q) [cm-1]
+        """
+        self._create_modelObject()
                 
         norm =  1.0e8/self.params['lores_density']*self.params['scale']
-        #return norm*pointsmodelpy.get_lores_i(self.lores_model, q)
         return norm*pointsmodelpy.get_complex_iq_2D(self.complex_model, self.points, qx, qy)\
             + self.params['background']
@@ -624 +633 @@
             self.getPr()
 
-        # By dividing by the density instead of the actuall V/N, 
+        # By dividing by the density instead of the actual V/N, 
         # we have an uncertainty of +-1 on N because the number 
         # of points chosen for the simulation is int(density*volume).
@@ -653 +662 @@
         simerr = 2*val/self.npts
         return val, err+simerr
-        
+
+    def getIq2DError(self, qx, qy):
+        """
+            Return the simulated value along with its estimated
+            error for a given q-value
+            
+            Propagation of errors is used to evaluate the
+            uncertainty.
+            
+            @param q: q-value [float]
+            @return: mean, error [float, float]
+        """
+        self._create_modelObject()
+                
+        norm =  1.0e8/self.params['lores_density']*self.params['scale']
+        val = norm*pointsmodelpy.get_complex_iq_2D(self.complex_model, self.points, qx, qy)\
+            + self.params['background']
+        
+        # Simulation error (statistical)
+        norm =  1.0e8/self.params['lores_density']*self.params['scale'] \
+                * math.pow(self.npts/self.params['lores_density'], 1.0/3.0)/self.npts
+        err = norm*pointsmodelpy.get_complex_iq_2D_err(self.complex_model, self.points, qx, qy)
+        # Error on V/N
+        simerr = 2*val/self.npts
+        
+        # The error used for the position is over-simplified.
+        # The actual error was empirically found to be about
+        # an order of magnitude larger.
+        return val, 10.0*err+simerr
+        

sansrealspace/src/realspace/test/early_test.py

@@ -237 +237 @@
     print ana, sim, sim/ana, ana/sim
     
-
+def test_7():
+    from sans.models.CoreShellModel import CoreShellModel
+    print "Testing core-shell"
+    radius = 15
+    thickness = 5
+    density = 5
+    
+    core_vol = 4.0/3.0*math.pi*radius*radius*radius
+    outer_radius = radius+thickness
+    shell_vol = 4.0/3.0*math.pi*outer_radius*outer_radius*outer_radius - core_vol
+    shell_sld = -1.0*core_vol/shell_vol
+
+    # Core-shell
+    sphere = CoreShellModel()
+    # Core radius
+    sphere.setParam('radius', radius)
+    # Shell thickness
+    sphere.setParam('thickness', thickness)
+    sphere.setParam('core_sld', 1.0)
+    sphere.setParam('shell_sld', shell_sld)
+    sphere.setParam('solvent_sld', 0.0)
+    sphere.setParam('background', 0.0)
+    sphere.setParam('scale', 1.0)
+    ana = sphere
+   
+    canvas = VolumeCanvas.VolumeCanvas()        
+    canvas.setParam('lores_density', density)
+    
+    handle = canvas.add('sphere')
+    canvas.setParam('%s.radius' % handle, outer_radius)
+    canvas.setParam('%s.contrast' % handle, shell_sld)
+   
+    handle2 = canvas.add('sphere')
+    canvas.setParam('%s.radius' % handle2, radius)
+    canvas.setParam('%s.contrast' % handle2, 1.0)
+           
+    canvas.setParam('scale' , 1.0)
+    canvas.setParam('background' , 0.0)
+    
+               
+    """ Testing default core-shell orientation """
+    qlist = [.0001, 0.002, .01, .1, 1.0, 5.]
+    for q in qlist:
+        ana_val = ana.runXY([q, 0.2])
+        sim_val, err = canvas.getIq2DError(q, 0.2)
+        print ana_val, sim_val, sim_val/ana_val, err, (sim_val-ana_val)/err
+    
+
+    
     
 if __name__ == "__main__":
-    test_6()
+    test_7()

sansrealspace/src/realspace/test/sim_validation.py

@@ -18 +18 @@
     def __init__(self):
         self.density = 0.1
+        self.canvas  = None
+        self.ana     = None
         self.create()
         
@@ -23 +25 @@
         pass
          
-    def run_sim(self, q, density=None):
+    def run_sim2D(self, qx, qy, density=None):
         """
             Calculate the mean and error of the simulation
@@ -34 +36 @@
             self.create()
         
+        return self.canvas.getIq2DError(qx, qy)
+    
+    def run_sim(self, q, density=None):
+        """
+            Calculate the mean and error of the simulation
+            @param q: q-value to calculate at
+            @param density: point density of simulation
+            #return: mean, error
+        """
+        if not density == None:
+            self.density = density
+            self.create()
+        
         return self.canvas.getIqError(q)
+         
+    def run_ana2D(self, qx, qy):
+        """
+            Return analytical value
+            @param q: q-value to evaluate at [float]
+            @return: analytical output [float]
+        """
+        return self.ana.runXY([qx, qy]) 
          
     def run_ana(self, q):
@@ -227 +250 @@
     output.close()
        
+def validate_model_2D(validator, q_min, q_max, phi, n_q):
+    """
+         Validate a model
+         An output file containing a comparison between
+         simulation and the analytical solution will be
+         produced.
+         
+         @param validator: validator object
+         @param q_min: minimum q
+         @param q_max: maximum q
+         @param n_q: number of q points
+         @param N: number of times to evaluate each simulation point
+    """
+    
+    q_list = pylab.arange(q_min, q_max*1.0001, (q_max-q_min)/(n_q-1))
+    
+    output = open('%s_d=%g_Iq2D.txt' % (validator.name, validator.density), 'w')
+    output.write("PARS: %s\n" % validator.ana.params)
+    output.write("<q>  <ana>  <sim>  <err>\n")
+    t_0 = time.time()
+    for q in q_list:
+        ana = validator.run_ana2D(q*math.cos(phi), q*math.sin(phi))
+        sim, err = validator.run_sim2D(q*math.cos(phi), q*math.sin(phi))
+        print "q=%-g  ana=%-g  sim=%-g  err=%-g  diff=%-g (%-g) %s" % (q, ana, sim, err, 
+                        (sim-ana), sim/ana, str(math.fabs(sim-ana)>err))
+        output.write("%g  %g  %g  %g\n" % (q, ana, sim, err))
+    print "Time elapsed: ", time.time()-t_0
+    output.close()
+       
 def check_density(validator, q, d_min, d_max, n_d): 
     """
@@ -257 +309 @@
     
 if __name__ == '__main__':
-    #vali = CoreShellValidator(radius = 15, thickness=5, density = 0.01)
-    #validate_model(vali, q_min=0.001, q_max=1, n_q=50, N=40)
-
+    
+    # 2D: Density=5, 71.2 secs for 50 points
+    vali = CoreShellValidator(radius = 15, thickness=5, density = 5.0)
+    #validate_model(vali, q_min=0.001, q_max=1, n_q=50)
+    validate_model_2D(vali, q_min=0.001, q_max=1, phi=1.0, n_q=50)
+
+    # 2D: Density=2, 11.1 secs for 25 points
     #vali = SphereValidator(radius = 20, density = 0.02)
-    #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25, N=20)
-    
-    vali = CylinderValidator(radius = 20, length=100, density = 0.1)
-    validate_model(vali, q_min=0.001, q_max=0.5, n_q=25)
-    
+    #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25)
+    #vali = SphereValidator(radius = 20, density = 2.0)
+    #validate_model_2D(vali, q_min=0.001, q_max=0.5, phi=1.0, n_q=25)
+    
+    # 2D: Density=1, 19.4 secs for 25 points
+    # 2D: Density=0.5, 9.8 secs for 25 points
+    #vali = CylinderValidator(radius = 20, length=100, density = 0.1)
+    #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25)
+    #vali = CylinderValidator(radius = 20, length=100, density = 0.5)
+    #validate_model_2D(vali, q_min=0.001, q_max=0.2, phi=1.0, n_q=25)
+    
+    # 2D: Density=0.5, 2.26 secs for 25 points
     #vali = EllipsoidValidator(radius_a = 20, radius_b=15, density = 0.05)
     #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25)
+    #vali = EllipsoidValidator(radius_a = 20, radius_b=15, density = 0.5)
+    #validate_model_2D(vali, q_min=0.001, q_max=0.5, phi=1.0, n_q=25)
     
     #vali = HelixValidator(density = 0.05)

sansrealspace/src/realspace/test/utest_oriented.py

@@ -240 +240 @@
             self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
         except:
-            print ana_val, sim_val, sim_val/ana_val
+            print "Error", ana_val, sim_val, sim_val/ana_val
             raise sys.exc_type, sys.exc_value
 
@@ -288 +288 @@
         canvas.setParam('background' , 0.0)
         self.canvas = canvas 
+           
+    def testdefault(self):
+        """ Testing default core-shell orientation """
+        ana_val = self.ana.runXY([0.1, 0.2])
+        sim_val, err = self.canvas.getIq2DError(0.1, 0.2)
+        
+        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+                   
+class TestCoreShellError(unittest.TestCase):
+    """ Tests for oriented (2D) systems """
+        
+    def setUp(self):
+        """ Set up zero-SLD-average core-shell model """
+        from sans.models.CoreShellModel import CoreShellModel
+        
+        radius = 15
+        thickness = 5
+        density = 5
+        
+        core_vol = 4.0/3.0*math.pi*radius*radius*radius
+        self.outer_radius = radius+thickness
+        shell_vol = 4.0/3.0*math.pi*self.outer_radius*self.outer_radius*self.outer_radius - core_vol
+        self.shell_sld = -1.0*core_vol/shell_vol
+
+        self.density = density
+           
+        # Core-shell
+        sphere = CoreShellModel()
+        # Core radius
+        sphere.setParam('radius', radius)
+        # Shell thickness
+        sphere.setParam('thickness', thickness)
+        sphere.setParam('core_sld', 1.0)
+        sphere.setParam('shell_sld', self.shell_sld)
+        sphere.setParam('solvent_sld', 0.0)
+        sphere.setParam('background', 0.0)
+        sphere.setParam('scale', 1.0)
+        self.ana = sphere
+       
+        canvas = VolumeCanvas.VolumeCanvas()        
+        canvas.setParam('lores_density', self.density)
+        
+        handle = canvas.add('sphere')
+        canvas.setParam('%s.radius' % handle, self.outer_radius)
+        canvas.setParam('%s.contrast' % handle, self.shell_sld)
+       
+        handle2 = canvas.add('sphere')
+        canvas.setParam('%s.radius' % handle2, radius)
+        canvas.setParam('%s.contrast' % handle2, 1.0)
+               
+        canvas.setParam('scale' , 1.0)
+        canvas.setParam('background' , 0.0)
+        self.canvas = canvas 
                    
     def testdefault(self):
         """ Testing default core-shell orientation """
         ana_val = self.ana.runXY([0.1, 0.2])
-        sim_val = self.canvas.getIq2D(0.1, 0.2)
-        #print ana_val, sim_val, sim_val/ana_val
-        
-        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        sim_val, err = self.canvas.getIq2DError(0.1, 0.2)
+        
+        self.assert_( math.fabs(sim_val-ana_val) < 3.0 * err )
 
 class TestRunMethods(unittest.TestCase):
@@ -306 +358 @@
         radius_a = 10
         radius_b = 15
-        density = 1
+        density = 5
         
         self.ana = EllipsoidModel()
@@ -339 +391 @@
         #print ana_val, sim_val, sim_val/ana_val
         
-        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        try:
+            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        except:
+            print "Error", ana_val, sim_val, sim_val/ana_val
+            raise sys.exc_type, sys.exc_value
 
     def testRunXY_float(self):
@@ -347 +403 @@
         #print ana_val, sim_val, sim_val/ana_val
         
-        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        try:
+            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        except:
+            print "Error", ana_val, sim_val, sim_val/ana_val
+            raise sys.exc_type, sys.exc_value
 
     def testRun_float(self):
@@ -355 +415 @@
         #print ana_val, sim_val, sim_val/ana_val
         
-        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        try:
+            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        except:
+            print "Error", ana_val, sim_val, sim_val/ana_val
+            raise sys.exc_type, sys.exc_value
 
     def testRun_list(self):
@@ -363 +427 @@
         #print ana_val, sim_val, sim_val/ana_val
         
-        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        try:
+            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 )
+        except:
+            print "Error", ana_val, sim_val, sim_val/ana_val
+            raise sys.exc_type, sys.exc_value