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