source: sasview/src/sas/sascalc/realspace/VolumeCanvas.py

ESS_GUI
Last change on this file was 98e3f24, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

code cleanup

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