Changeset a2c1196 in sasview for sansrealspace


Ignore:
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
Files:
10 added
5 edited

Legend:

Unmodified
Added
Removed
  • sansrealspace/.pydevproject

    rba1d1e9 ra2c1196  
    66<pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH"> 
    77<path>/RealSpaceTopLayer/src</path> 
     8<path>/RealSpaceTopLayer/src/realspace</path> 
    89</pydev_pathproperty> 
    910</pydev_project> 
  • sansrealspace/src/realspace/VolumeCanvas.py

    r3c75696 ra2c1196  
    551551            raise ValueError, "runXY(q): bad type for q" 
    552552     
    553          
    554     def getIq2D(self, qx, qy): 
     553    def _create_modelObject(self): 
    555554        """ 
    556555            Returns simulate I(q) for given q_x and q_y values. 
     556            Also returns model object 
    557557            @param qx: q_x [A-1] 
    558558            @param qy: q_y [A-1] 
    559             @return: I(q) [cm-1] 
     559            @return: I(q) [cm-1], model object 
    560560        """ 
    561561        # To find a complete example of the correct call order: 
     
    581581        #pointsmodelpy.get_lorespoints(self.lores_model, self.points) 
    582582        self.npts = pointsmodelpy.get_complexpoints(self.complex_model, self.points) 
     583         
     584         
     585    def getIq2D(self, qx, qy): 
     586        """ 
     587            Returns simulate I(q) for given q_x and q_y values. 
     588            @param qx: q_x [A-1] 
     589            @param qy: q_y [A-1] 
     590            @return: I(q) [cm-1] 
     591        """ 
     592        self._create_modelObject() 
    583593                 
    584594        norm =  1.0e8/self.params['lores_density']*self.params['scale'] 
    585         #return norm*pointsmodelpy.get_lores_i(self.lores_model, q) 
    586595        return norm*pointsmodelpy.get_complex_iq_2D(self.complex_model, self.points, qx, qy)\ 
    587596            + self.params['background'] 
     
    624633            self.getPr() 
    625634 
    626         # By dividing by the density instead of the actuall V/N,  
     635        # By dividing by the density instead of the actual V/N,  
    627636        # we have an uncertainty of +-1 on N because the number  
    628637        # of points chosen for the simulation is int(density*volume). 
     
    653662        simerr = 2*val/self.npts 
    654663        return val, err+simerr 
    655          
     664 
     665    def getIq2DError(self, qx, qy): 
     666        """ 
     667            Return the simulated value along with its estimated 
     668            error for a given q-value 
     669             
     670            Propagation of errors is used to evaluate the 
     671            uncertainty. 
     672             
     673            @param q: q-value [float] 
     674            @return: mean, error [float, float] 
     675        """ 
     676        self._create_modelObject() 
     677                 
     678        norm =  1.0e8/self.params['lores_density']*self.params['scale'] 
     679        val = norm*pointsmodelpy.get_complex_iq_2D(self.complex_model, self.points, qx, qy)\ 
     680            + self.params['background'] 
     681         
     682        # Simulation error (statistical) 
     683        norm =  1.0e8/self.params['lores_density']*self.params['scale'] \ 
     684                * math.pow(self.npts/self.params['lores_density'], 1.0/3.0)/self.npts 
     685        err = norm*pointsmodelpy.get_complex_iq_2D_err(self.complex_model, self.points, qx, qy) 
     686        # Error on V/N 
     687        simerr = 2*val/self.npts 
     688         
     689        # The error used for the position is over-simplified. 
     690        # The actual error was empirically found to be about 
     691        # an order of magnitude larger. 
     692        return val, 10.0*err+simerr 
     693         
  • sansrealspace/src/realspace/test/early_test.py

    r1b0707e1 ra2c1196  
    237237    print ana, sim, sim/ana, ana/sim 
    238238     
    239  
     239def test_7(): 
     240    from sans.models.CoreShellModel import CoreShellModel 
     241    print "Testing core-shell" 
     242    radius = 15 
     243    thickness = 5 
     244    density = 5 
     245     
     246    core_vol = 4.0/3.0*math.pi*radius*radius*radius 
     247    outer_radius = radius+thickness 
     248    shell_vol = 4.0/3.0*math.pi*outer_radius*outer_radius*outer_radius - core_vol 
     249    shell_sld = -1.0*core_vol/shell_vol 
     250 
     251    # Core-shell 
     252    sphere = CoreShellModel() 
     253    # Core radius 
     254    sphere.setParam('radius', radius) 
     255    # Shell thickness 
     256    sphere.setParam('thickness', thickness) 
     257    sphere.setParam('core_sld', 1.0) 
     258    sphere.setParam('shell_sld', shell_sld) 
     259    sphere.setParam('solvent_sld', 0.0) 
     260    sphere.setParam('background', 0.0) 
     261    sphere.setParam('scale', 1.0) 
     262    ana = sphere 
     263    
     264    canvas = VolumeCanvas.VolumeCanvas()         
     265    canvas.setParam('lores_density', density) 
     266     
     267    handle = canvas.add('sphere') 
     268    canvas.setParam('%s.radius' % handle, outer_radius) 
     269    canvas.setParam('%s.contrast' % handle, shell_sld) 
     270    
     271    handle2 = canvas.add('sphere') 
     272    canvas.setParam('%s.radius' % handle2, radius) 
     273    canvas.setParam('%s.contrast' % handle2, 1.0) 
     274            
     275    canvas.setParam('scale' , 1.0) 
     276    canvas.setParam('background' , 0.0) 
     277     
     278                
     279    """ Testing default core-shell orientation """ 
     280    qlist = [.0001, 0.002, .01, .1, 1.0, 5.] 
     281    for q in qlist: 
     282        ana_val = ana.runXY([q, 0.2]) 
     283        sim_val, err = canvas.getIq2DError(q, 0.2) 
     284        print ana_val, sim_val, sim_val/ana_val, err, (sim_val-ana_val)/err 
     285     
     286 
     287     
    240288     
    241289if __name__ == "__main__": 
    242     test_6() 
     290    test_7() 
  • sansrealspace/src/realspace/test/sim_validation.py

    rba1d1e9 ra2c1196  
    1818    def __init__(self): 
    1919        self.density = 0.1 
     20        self.canvas  = None 
     21        self.ana     = None 
    2022        self.create() 
    2123         
     
    2325        pass 
    2426          
    25     def run_sim(self, q, density=None): 
     27    def run_sim2D(self, qx, qy, density=None): 
    2628        """ 
    2729            Calculate the mean and error of the simulation 
     
    3436            self.create() 
    3537         
     38        return self.canvas.getIq2DError(qx, qy) 
     39     
     40    def run_sim(self, q, density=None): 
     41        """ 
     42            Calculate the mean and error of the simulation 
     43            @param q: q-value to calculate at 
     44            @param density: point density of simulation 
     45            #return: mean, error 
     46        """ 
     47        if not density == None: 
     48            self.density = density 
     49            self.create() 
     50         
    3651        return self.canvas.getIqError(q) 
     52          
     53    def run_ana2D(self, qx, qy): 
     54        """ 
     55            Return analytical value 
     56            @param q: q-value to evaluate at [float] 
     57            @return: analytical output [float] 
     58        """ 
     59        return self.ana.runXY([qx, qy])  
    3760          
    3861    def run_ana(self, q): 
     
    227250    output.close() 
    228251        
     252def validate_model_2D(validator, q_min, q_max, phi, n_q): 
     253    """ 
     254         Validate a model 
     255         An output file containing a comparison between 
     256         simulation and the analytical solution will be 
     257         produced. 
     258          
     259         @param validator: validator object 
     260         @param q_min: minimum q 
     261         @param q_max: maximum q 
     262         @param n_q: number of q points 
     263         @param N: number of times to evaluate each simulation point 
     264    """ 
     265     
     266    q_list = pylab.arange(q_min, q_max*1.0001, (q_max-q_min)/(n_q-1)) 
     267     
     268    output = open('%s_d=%g_Iq2D.txt' % (validator.name, validator.density), 'w') 
     269    output.write("PARS: %s\n" % validator.ana.params) 
     270    output.write("<q>  <ana>  <sim>  <err>\n") 
     271    t_0 = time.time() 
     272    for q in q_list: 
     273        ana = validator.run_ana2D(q*math.cos(phi), q*math.sin(phi)) 
     274        sim, err = validator.run_sim2D(q*math.cos(phi), q*math.sin(phi)) 
     275        print "q=%-g  ana=%-g  sim=%-g  err=%-g  diff=%-g (%-g) %s" % (q, ana, sim, err,  
     276                        (sim-ana), sim/ana, str(math.fabs(sim-ana)>err)) 
     277        output.write("%g  %g  %g  %g\n" % (q, ana, sim, err)) 
     278    print "Time elapsed: ", time.time()-t_0 
     279    output.close() 
     280        
    229281def check_density(validator, q, d_min, d_max, n_d):  
    230282    """ 
     
    257309     
    258310if __name__ == '__main__': 
    259     #vali = CoreShellValidator(radius = 15, thickness=5, density = 0.01) 
    260     #validate_model(vali, q_min=0.001, q_max=1, n_q=50, N=40) 
    261  
     311     
     312    # 2D: Density=5, 71.2 secs for 50 points 
     313    vali = CoreShellValidator(radius = 15, thickness=5, density = 5.0) 
     314    #validate_model(vali, q_min=0.001, q_max=1, n_q=50) 
     315    validate_model_2D(vali, q_min=0.001, q_max=1, phi=1.0, n_q=50) 
     316 
     317    # 2D: Density=2, 11.1 secs for 25 points 
    262318    #vali = SphereValidator(radius = 20, density = 0.02) 
    263     #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25, N=20) 
    264      
    265     vali = CylinderValidator(radius = 20, length=100, density = 0.1) 
    266     validate_model(vali, q_min=0.001, q_max=0.5, n_q=25) 
    267      
     319    #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25) 
     320    #vali = SphereValidator(radius = 20, density = 2.0) 
     321    #validate_model_2D(vali, q_min=0.001, q_max=0.5, phi=1.0, n_q=25) 
     322     
     323    # 2D: Density=1, 19.4 secs for 25 points 
     324    # 2D: Density=0.5, 9.8 secs for 25 points 
     325    #vali = CylinderValidator(radius = 20, length=100, density = 0.1) 
     326    #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25) 
     327    #vali = CylinderValidator(radius = 20, length=100, density = 0.5) 
     328    #validate_model_2D(vali, q_min=0.001, q_max=0.2, phi=1.0, n_q=25) 
     329     
     330    # 2D: Density=0.5, 2.26 secs for 25 points 
    268331    #vali = EllipsoidValidator(radius_a = 20, radius_b=15, density = 0.05) 
    269332    #validate_model(vali, q_min=0.001, q_max=0.5, n_q=25) 
     333    #vali = EllipsoidValidator(radius_a = 20, radius_b=15, density = 0.5) 
     334    #validate_model_2D(vali, q_min=0.001, q_max=0.5, phi=1.0, n_q=25) 
    270335     
    271336    #vali = HelixValidator(density = 0.05) 
  • sansrealspace/src/realspace/test/utest_oriented.py

    r3c75696 ra2c1196  
    240240            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
    241241        except: 
    242             print ana_val, sim_val, sim_val/ana_val 
     242            print "Error", ana_val, sim_val, sim_val/ana_val 
    243243            raise sys.exc_type, sys.exc_value 
    244244 
     
    288288        canvas.setParam('background' , 0.0) 
    289289        self.canvas = canvas  
     290            
     291    def testdefault(self): 
     292        """ Testing default core-shell orientation """ 
     293        ana_val = self.ana.runXY([0.1, 0.2]) 
     294        sim_val, err = self.canvas.getIq2DError(0.1, 0.2) 
     295         
     296        self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     297                    
     298class TestCoreShellError(unittest.TestCase): 
     299    """ Tests for oriented (2D) systems """ 
     300         
     301    def setUp(self): 
     302        """ Set up zero-SLD-average core-shell model """ 
     303        from sans.models.CoreShellModel import CoreShellModel 
     304         
     305        radius = 15 
     306        thickness = 5 
     307        density = 5 
     308         
     309        core_vol = 4.0/3.0*math.pi*radius*radius*radius 
     310        self.outer_radius = radius+thickness 
     311        shell_vol = 4.0/3.0*math.pi*self.outer_radius*self.outer_radius*self.outer_radius - core_vol 
     312        self.shell_sld = -1.0*core_vol/shell_vol 
     313 
     314        self.density = density 
     315            
     316        # Core-shell 
     317        sphere = CoreShellModel() 
     318        # Core radius 
     319        sphere.setParam('radius', radius) 
     320        # Shell thickness 
     321        sphere.setParam('thickness', thickness) 
     322        sphere.setParam('core_sld', 1.0) 
     323        sphere.setParam('shell_sld', self.shell_sld) 
     324        sphere.setParam('solvent_sld', 0.0) 
     325        sphere.setParam('background', 0.0) 
     326        sphere.setParam('scale', 1.0) 
     327        self.ana = sphere 
     328        
     329        canvas = VolumeCanvas.VolumeCanvas()         
     330        canvas.setParam('lores_density', self.density) 
     331         
     332        handle = canvas.add('sphere') 
     333        canvas.setParam('%s.radius' % handle, self.outer_radius) 
     334        canvas.setParam('%s.contrast' % handle, self.shell_sld) 
     335        
     336        handle2 = canvas.add('sphere') 
     337        canvas.setParam('%s.radius' % handle2, radius) 
     338        canvas.setParam('%s.contrast' % handle2, 1.0) 
     339                
     340        canvas.setParam('scale' , 1.0) 
     341        canvas.setParam('background' , 0.0) 
     342        self.canvas = canvas  
    290343                    
    291344    def testdefault(self): 
    292345        """ Testing default core-shell orientation """ 
    293346        ana_val = self.ana.runXY([0.1, 0.2]) 
    294         sim_val = self.canvas.getIq2D(0.1, 0.2) 
    295         #print ana_val, sim_val, sim_val/ana_val 
    296          
    297         self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     347        sim_val, err = self.canvas.getIq2DError(0.1, 0.2) 
     348         
     349        self.assert_( math.fabs(sim_val-ana_val) < 3.0 * err ) 
    298350 
    299351class TestRunMethods(unittest.TestCase): 
     
    306358        radius_a = 10 
    307359        radius_b = 15 
    308         density = 1 
     360        density = 5 
    309361         
    310362        self.ana = EllipsoidModel() 
     
    339391        #print ana_val, sim_val, sim_val/ana_val 
    340392         
    341         self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     393        try: 
     394            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     395        except: 
     396            print "Error", ana_val, sim_val, sim_val/ana_val 
     397            raise sys.exc_type, sys.exc_value 
    342398 
    343399    def testRunXY_float(self): 
     
    347403        #print ana_val, sim_val, sim_val/ana_val 
    348404         
    349         self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     405        try: 
     406            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     407        except: 
     408            print "Error", ana_val, sim_val, sim_val/ana_val 
     409            raise sys.exc_type, sys.exc_value 
    350410 
    351411    def testRun_float(self): 
     
    355415        #print ana_val, sim_val, sim_val/ana_val 
    356416         
    357         self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     417        try: 
     418            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     419        except: 
     420            print "Error", ana_val, sim_val, sim_val/ana_val 
     421            raise sys.exc_type, sys.exc_value 
    358422 
    359423    def testRun_list(self): 
     
    363427        #print ana_val, sim_val, sim_val/ana_val 
    364428         
    365         self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     429        try: 
     430            self.assert_( math.fabs(sim_val/ana_val-1.0)<0.05 ) 
     431        except: 
     432            print "Error", ana_val, sim_val, sim_val/ana_val 
     433            raise sys.exc_type, sys.exc_value 
    366434 
    367435           
Note: See TracChangeset for help on using the changeset viewer.