source: sasview/sansrealspace/src/realspace/VolumeCanvas.py @ a2c1196

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 a2c1196 was a2c1196, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

Added 2D (oriented systems) error estimation, tests and validation.
  • Property mode set to 100644
File size: 23.7 KB
Line 
1#!/usr/bin/env python
2""" Volume Canvas
3    Simulation canvas
4"""
5
6from sans.models.BaseComponent import BaseComponent
7from sansModeling.pointsmodelpy import pointsmodelpy
8from sansModeling.geoshapespy import geoshapespy
9
10import os.path, math
11
12class ShapeDescriptor:
13    """
14        Class to hold the information about a shape
15    """
16    def __init__(self):
17        """
18            Initialization
19        """
20        ## Real space object
21        self.shapeObject = None
22        ## Parameters of the object
23        self.params = {}
24        self.params["center"] = [0, 0, 0]
25        # Orientation are angular offsets in degrees with respect to X, Y, Z
26        self.params["orientation"] = [0, 0, 0]
27        # Default to lores shape
28        self.params['is_lores'] = True
29        self.params['order'] = 0
30           
31    def create(self):
32        """
33            Create an instance of the shape
34        """
35        # Set center
36        x0 = self.params["center"][0]
37        y0 = self.params["center"][1]
38        z0 = self.params["center"][2]
39        geoshapespy.set_center(self.shapeObject, x0, y0, z0)
40       
41        # Set orientation
42        x0 = self.params["orientation"][0]
43        y0 = self.params["orientation"][1]
44        z0 = self.params["orientation"][2]
45        geoshapespy.set_orientation(self.shapeObject, x0, y0, z0)
46               
47class SphereDescriptor(ShapeDescriptor):
48    """
49        Descriptor for a sphere
50    """
51    def __init__(self):
52        """
53            Initialization
54        """ 
55        ShapeDescriptor.__init__(self)
56        # Default parameters
57        self.params["type"]   = "sphere"
58        self.params["radius"] = 20.0
59        self.params["contrast"] = 1.0
60
61    def create(self):
62        """
63            Create an instance of the shape
64            @return: instance of the shape
65        """
66        self.shapeObject = geoshapespy.new_sphere(\
67            self.params["radius"])
68       
69        ShapeDescriptor.create(self)   
70        return self.shapeObject
71   
72class CylinderDescriptor(ShapeDescriptor):
73    """
74        Descriptor for a cylinder
75        Orientation: Default cylinder is along Y
76    """
77    def __init__(self):
78        """
79            Initialization
80        """ 
81        ShapeDescriptor.__init__(self)
82        # Default parameters
83        self.params["type"]   = "cylinder"
84        self.params["length"] = 40.0
85        self.params["radius"] = 10.0
86        self.params["contrast"] = 1.0
87       
88    def create(self):
89        """
90            Create an instance of the shape
91            @return: instance of the shape
92        """
93        self.shapeObject = geoshapespy.new_cylinder(\
94            self.params["radius"], self.params["length"])
95
96        ShapeDescriptor.create(self)
97        return self.shapeObject
98       
99
100class EllipsoidDescriptor(ShapeDescriptor):
101    """
102        Descriptor for an ellipsoid
103    """
104    def __init__(self):
105        """
106            Initialization
107        """ 
108        ShapeDescriptor.__init__(self)
109        # Default parameters
110        self.params["type"]   = "ellipsoid"
111        self.params["radius_x"] = 30.0
112        self.params["radius_y"] = 20.0
113        self.params["radius_z"] = 10.0
114        self.params["contrast"] = 1.0
115       
116    def create(self):
117        """
118            Create an instance of the shape
119            @return: instance of the shape
120        """
121        self.shapeObject = geoshapespy.new_ellipsoid(\
122            self.params["radius_x"], self.params["radius_y"], 
123            self.params["radius_z"])
124       
125        ShapeDescriptor.create(self)   
126        return self.shapeObject
127       
128class HelixDescriptor(ShapeDescriptor):
129    """
130        Descriptor for an helix
131    """
132    def __init__(self):
133        """
134            Initialization
135        """ 
136        ShapeDescriptor.__init__(self)
137        # Default parameters
138        self.params["type"]   = "singlehelix"
139        self.params["radius_helix"] = 10.0
140        self.params["radius_tube"] = 3.0
141        self.params["pitch"] = 34.0
142        self.params["turns"] = 3.0
143        self.params["contrast"] = 1.0
144
145    def create(self):
146        """
147            Create an instance of the shape
148            @return: instance of the shape
149        """
150        self.shapeObject = geoshapespy.new_singlehelix(\
151            self.params["radius_helix"], self.params["radius_tube"], 
152            self.params["pitch"], self.params["turns"])
153       
154        ShapeDescriptor.create(self)   
155        return self.shapeObject
156       
157class PDBDescriptor(ShapeDescriptor):
158    """
159        Descriptor for a PDB set of points
160    """
161    def __init__(self, filename):
162        """
163            Initialization
164            @param filename: name of the PDB file to load
165        """ 
166        ShapeDescriptor.__init__(self)
167        # Default parameters
168        self.params["type"]   = "pdb"
169        self.params["file"] = filename
170        self.params['is_lores'] = False
171
172    def create(self):
173        """
174            Create an instance of the shape
175            @return: instance of the shape
176        """
177        self.shapeObject = pointsmodelpy.new_pdbmodel()
178        pointsmodelpy.pdbmodel_add(self.shapeObject, self.params['file'])       
179       
180        #ShapeDescriptor.create(self)   
181        return self.shapeObject
182       
183# Define a dictionary for the shape until we find
184# a better way to create them
185shape_dict = {'sphere':SphereDescriptor,
186              'cylinder':CylinderDescriptor,
187              'ellipsoid':EllipsoidDescriptor,
188              'singlehelix':HelixDescriptor}
189       
190class VolumeCanvas(BaseComponent):
191    """
192        Class representing an empty space volume to add
193        geometrical object to.
194    """
195   
196    def __init__(self):
197        """
198            Initialization
199        """
200        BaseComponent.__init__(self)
201       
202        ## Maximum value of q reachable
203        self.params['q_max'] = 0.1
204        self.params['lores_density'] = 0.1
205        self.params['scale'] = 1.0
206        self.params['background'] = 0.0
207       
208        self.lores_model = pointsmodelpy.new_loresmodel(self.params['lores_density'])
209        self.complex_model = pointsmodelpy.new_complexmodel()
210        self.shapes = {}
211        self.shapecount = 0       
212        self.points = None
213        self.npts = 0
214        self.hasPr = False       
215       
216    def addObject(self, shapeDesc, id = None):
217        """
218            Adds a real-space object to the canvas.
219       
220            @param shapeDesc: object to add to the canvas [ShapeDescriptor]
221            @param id: string handle for the object [string] [optional]
222            @return: string handle for the object
223        """
224        # If the handle is not provided, create one
225        if id == None:
226            id = shapeDesc.params["type"]+str(self.shapecount)
227         
228        # Self the order number
229        shapeDesc.params['order'] = self.shapecount
230        # Store the shape in a dictionary entry associated
231        # with the handle
232        self.shapes[id] = shapeDesc
233        self.shapecount += 1
234
235        #model changed, need to recalculate P(r)
236        self.hasPr = False
237
238        return id
239           
240   
241    def add(self, shape, id = None):
242        """
243            The intend of this method is to eventually be able to use it
244            as a factory for the canvas and unify the simulation with the
245            analytical solutions. For instance, if one adds a cylinder and
246            it is the only shape on the canvas, the analytical solution
247            could be called. If multiple shapes are involved, then
248            simulation has to be performed.
249       
250            @param shape: name of the object to add to the canvas [string]
251            @param id: string handle for the object [string] [optional]
252            @return: string handle for the object
253        """
254        # If the handle is not provided, create one
255        if id == None:
256            id = "shape"+str(self.shapecount)
257 
258        #shapeDesc = ShapeDescriptor(shape.lower())
259        if shape.lower() in shape_dict:
260            shapeDesc = shape_dict[shape.lower()]()
261        elif os.path.isfile(shape):
262            # A valid filename was supplier, create a PDB object
263            shapeDesc = PDBDescriptor(shape)
264        else:
265            raise ValueError, "VolumeCanvas.add: Unknown shape %s" % shape
266       
267        # Self the order number
268        shapeDesc.params['order'] = self.shapecount
269        # Store the shape in a dictionary entry associated
270        # with the handle
271        self.shapes[id] = shapeDesc
272        self.shapecount += 1
273
274        #model changed, need to recalculate P(r)
275        self.hasPr = False
276
277        return id
278
279    def delete(self, id):
280        """
281            Delete a shape. The ID for the shape is required.
282            @param id: string handle for the object [string] [optional]
283        """
284
285        if self.shapes.has_key(id):
286            del self.shapes[id]
287        else:
288            raise KeyError, "VolumeCanvas.delete: could not find shape ID"
289
290        #model changed, need to recalculate P(r)
291        self.hasPr = False
292
293
294    def setParam(self, name, value):   
295        """
296            @param name: name of the parameter to change
297            @param value: value to give the parameter
298        """
299       
300        # Lowercase for case insensitivity
301        name = name.lower()
302       
303        # Look for shape access
304        toks = name.split('.')
305       
306        # If a shape identifier was given, look the shape up
307        # in the dictionary
308        if len(toks)>1:
309            if toks[0] in self.shapes.keys():
310                # The shape was found, now look for the parameter
311                if toks[1] in self.shapes[toks[0]].params:
312                    # The parameter was found, now change it
313                    self.shapes[toks[0]].params[toks[1]] = value
314                    self.hasPr = False
315                else:
316                    raise ValueError, "Could not find parameter %s" % name
317            else:
318                raise ValueError, "Could not find shape %s" % toks[0]
319       
320        else:
321            # If we are not accessing the parameters of a
322            # shape, see if the parameter is part of this object
323            BaseComponent.setParam(self, name, value)
324            self.hasPr = False
325
326    def getParam(self, name):   
327        """
328            @param name: name of the parameter to change
329        """
330        #TODO: clean this up
331       
332        # Lowercase for case insensitivity
333        name = name.lower()
334       
335        # Look for sub-model access
336        toks = name.split('.')
337        if len(toks) == 1:
338            try:
339                self.params.has_key(toks[0])
340            except KeyError:
341                raise ValueError, \
342                    "VolumeCanvas.getParam: Could not find %s" % name
343
344            value = self.params[toks[0]]
345            if isinstance(value, ShapeDescriptor):
346                raise ValueError, \
347                    "VolumeCanvas.getParam: Cannot get parameter value." 
348            else:
349                return value
350
351        elif len(toks) == 2:
352            try:
353                self.shapes.has_key(toks[0])
354            except KeyError:
355                raise ValueError, \
356                    "VolumeCanvas.getParam: Could not find %s" % name
357
358            shapeinstance = self.shapes[toks[0]]
359
360            try:
361                shapeinstance.params.has_key(toks[1])
362            except KeyError:
363                raise ValueError, \
364                    "VolumeCanvas.getParam: Could not find %s" % name
365
366            return shapeinstance.params[toks[1]]
367
368        else:
369            raise ValueError, \
370                "VolumeCanvas.getParam: Could not find %s" % name
371           
372    def getParamList(self, shapeid = None):
373        """
374               return a full list of all available parameters from
375           self.params.keys(). If a key in self.params is a instance
376           of ShapeDescriptor, extend the return list to:
377           [param1,param2,shapeid.param1,shapeid.param2.......]
378
379           If shapeid is provided, return the list of parameters that
380           belongs to that shape id only : [shapeid.param1, shapeid.param2...]
381        """
382
383        param_list = []
384        if shapeid == None:       
385            for key1 in self.params.keys():
386                #value1 = self.params[key1]
387                param_list.append(key1)
388            for key2 in self.shapes.keys():
389                value2 = self.shapes[key2]
390                header = key2 + '.'
391                for key3 in value2.params.keys():   
392                    fullname = header + key3                 
393                    param_list.append(fullname)
394     
395        else:
396            try:
397                self.shapes.has_key(shapeid)
398            except KeyError:
399                raise ValueError, \
400                    "VolumeCanvas: getParamList: Could not find %s" % shapeid
401            header = shapeid + '.'
402            param_list = self.shapes[shapeid].params.keys() 
403            for i in range(len(param_list)):
404                param_list[i] = header + param_list[i]
405
406        return param_list
407
408    def getShapeList(self):
409        """
410            Return a list of the shapes
411        """
412        return self.shapes.keys()
413
414    def addSingleShape(self, shapeDesc):
415        """
416            create shapeobject based on shapeDesc
417            @param shapeDesc: shape description
418        """
419        #Create the object model
420        shapeDesc.create()
421                   
422        if shapeDesc.params['is_lores']:
423            # Add the shape to the lores_model
424            pointsmodelpy.lores_add(self.lores_model, 
425                shapeDesc.shapeObject, shapeDesc.params['contrast']) 
426
427    def createVolumeFromList(self):
428        """
429            Create a new lores model with all the shapes in our internal list
430            Whenever we change a parameter of a shape, we have to re-create
431            the whole thing.
432           
433            Items with higher 'order' number take precedence for regions
434            of space that are shared with other objects. Points in the
435            overlapping region belonging to objects with lower 'order'
436            will be ignored.
437           
438            Items are added in decreasing 'order' number.
439            The item with the highest 'order' will be added *first*.
440            [That conventions is prescribed by the realSpaceModeling module]
441        """
442       
443        # Create empty model
444        self.lores_model = \
445            pointsmodelpy.new_loresmodel(self.params['lores_density'])
446
447        # Create empty complex model
448        self.complex_model = pointsmodelpy.new_complexmodel()
449       
450        # Order the object first
451        obj_list = []
452   
453        for shape in self.shapes:
454            order = self.shapes[shape].params['order']
455            # find where to place it in the list
456            stored = False
457           
458            for i in range(len(obj_list)):
459                if obj_list[i][0] > order:
460                    obj_list.insert(i, [order, shape])
461                    stored = True
462                    break
463           
464            if not stored:
465                obj_list.append([order, shape])
466               
467        # Add each shape
468        len_list = len(obj_list)
469        for i in range(len_list-1, -1, -1):
470            shapedesc = self.shapes[obj_list[i][1]]
471            self.addSingleShape(shapedesc)
472
473        return 0     
474   
475    def getPr(self):
476        """
477            Calculate P(r)
478            @return: calculation output flag
479        """
480        # To find a complete example of the correct call order:
481        # In LORES2, in actionclass.py, method CalculateAction._get_iq()
482       
483        # If there are not shapes, do nothing
484        if len(self.shapes) == 0:
485            self.hasPr = False
486            return 0
487       
488        # generate space filling points from shape list
489        self.createVolumeFromList()
490
491        self.points = pointsmodelpy.new_point3dvec()
492
493        pointsmodelpy.complexmodel_add(self.complex_model, 
494                                        self.lores_model, "LORES")
495        for shape in self.shapes:
496            if self.shapes[shape].params['is_lores'] == False:
497                pointsmodelpy.complexmodel_add(self.complex_model, 
498                    self.shapes[shape].shapeObject, "PDB")
499       
500        #pointsmodelpy.get_lorespoints(self.lores_model, self.points)
501        self.npts = pointsmodelpy.get_complexpoints(self.complex_model, self.points)
502       
503        # expecting the rmax is a positive float or 0. The maximum distance.
504        #rmax = pointsmodelpy.get_lores_pr(self.lores_model, self.points)   
505         
506        rmax = pointsmodelpy.get_complex_pr(self.complex_model, self.points) 
507        self.hasPr = True   
508
509        return rmax
510       
511    def run(self, q = 0):
512        """
513            Returns the value of I(q) for a given q-value
514            @param q: q-value ([float] or [list]) ([A-1] or [[A-1], [rad]])
515            @return: I(q) [float] [cm-1]
516        """
517        # Check for 1D q length
518        if q.__class__.__name__ == 'int' \
519            or q.__class__.__name__ == 'float':
520            return self.getIq(q)
521        # Check for 2D q-value
522        elif q.__class__.__name__ == 'list':
523            # Compute (Qx, Qy) from (Q, phi)
524            # Phi is in radian and Q-values are in A-1
525            qx = q[0]*math.cos(q[1])
526            qy = q[0]*math.sin(q[1])
527            return self.getIq2D(qx, qy)
528        # Through an exception if it's not a
529        # type we recognize
530        else:
531            raise ValueError, "run(q): bad type for q"
532   
533    def runXY(self, q = 0):
534        """
535            Standard run command for the canvas.
536            Redirects to the correct method
537            according to the input type.
538            @param q: q-value [float] or [list] [A-1]
539            @return: I(q) [float] [cm-1]
540        """
541        # Check for 1D q length
542        if q.__class__.__name__ == 'int' \
543            or q.__class__.__name__ == 'float':
544            return self.getIq(q)
545        # Check for 2D q-value
546        elif q.__class__.__name__ == 'list':
547            return self.getIq2D(q[0], q[1])
548        # Through an exception if it's not a
549        # type we recognize
550        else:
551            raise ValueError, "runXY(q): bad type for q"
552   
553    def _create_modelObject(self):
554        """
555            Returns simulate I(q) for given q_x and q_y values.
556            Also returns model object
557            @param qx: q_x [A-1]
558            @param qy: q_y [A-1]
559            @return: I(q) [cm-1], model object
560        """
561        # To find a complete example of the correct call order:
562        # In LORES2, in actionclass.py, method CalculateAction._get_iq()
563       
564        # If there are not shapes, do nothing
565        if len(self.shapes) == 0:
566            self.hasPr = False
567            return 0
568       
569        # generate space filling points from shape list
570        self.createVolumeFromList()
571
572        self.points = pointsmodelpy.new_point3dvec()
573
574        pointsmodelpy.complexmodel_add(self.complex_model, 
575                                        self.lores_model, "LORES")
576        for shape in self.shapes:
577            if self.shapes[shape].params['is_lores'] == False:
578                pointsmodelpy.complexmodel_add(self.complex_model, 
579                    self.shapes[shape].shapeObject, "PDB")
580       
581        #pointsmodelpy.get_lorespoints(self.lores_model, self.points)
582        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()
593               
594        norm =  1.0e8/self.params['lores_density']*self.params['scale']
595        return norm*pointsmodelpy.get_complex_iq_2D(self.complex_model, self.points, qx, qy)\
596            + self.params['background']
597               
598       
599    def getIq(self, q):
600        """
601            Returns the value of I(q) for a given q-value
602           
603            This method should remain internal to the class
604            and the run() method should be used instead.
605           
606            @param q: q-value [float]
607            @return: I(q) [float]
608        """
609       
610        if self.hasPr == False:
611            self.getPr()
612
613        # By dividing by the density instead of the actuall V/N,
614        # we have an uncertainty of +-1 on N because the number
615        # of points chosen for the simulation is int(density*volume).
616        # Propagation of error gives:
617        #   delta(1/density^2) = 2*(1/density^2)/N
618        # where N is stored in self.npts
619
620        norm =  1.0e8/self.params['lores_density']*self.params['scale']
621        #return norm*pointsmodelpy.get_lores_i(self.lores_model, q)
622        return norm*pointsmodelpy.get_complex_i(self.complex_model, q)\
623            + self.params['background']
624   
625    def getError(self, q):
626        """
627            Returns the error of I(q) for a given q-value
628            @param q: q-value [float]
629            @return: I(q) [float]
630        """
631       
632        if self.hasPr == False:
633            self.getPr()
634
635        # By dividing by the density instead of the actual V/N,
636        # we have an uncertainty of +-1 on N because the number
637        # of points chosen for the simulation is int(density*volume).
638        # Propagation of error gives:
639        #   delta(1/density^2) = 2*(1/density^2)/N
640        # where N is stored in self.npts
641
642        norm =  1.0e8/self.params['lores_density']*self.params['scale']
643        #return norm*pointsmodelpy.get_lores_i(self.lores_model, q)
644        return norm*pointsmodelpy.get_complex_i_error(self.complex_model, q)\
645            + self.params['background']
646   
647    def getIqError(self, q):
648        """
649            Return the simulated value along with its estimated
650            error for a given q-value
651           
652            Propagation of errors is used to evaluate the
653            uncertainty.
654           
655            @param q: q-value [float]
656            @return: mean, error [float, float]
657        """
658        val = self.getIq(q)
659        # Simulation error (statistical)
660        err = self.getError(q)
661        # Error on V/N
662        simerr = 2*val/self.npts
663        return val, err+simerr
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       
Note: See TracBrowser for help on using the repository browser.