/[escript]/trunk/escript/py_src/modelframe.py
ViewVC logotype

Contents of /trunk/escript/py_src/modelframe.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2625 - (show annotations)
Fri Aug 21 06:30:25 2009 UTC (10 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 47318 byte(s)
Merging changes from new doco branch.
All docstrings are now in reStructured text.
A few email addresses have been fixes as well.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 Environment for implementing models in escript
24
25 :var __author__: name of author
26 :var __copyright__: copyrights
27 :var __license__: licence agreement
28 :var __url__: url entry point on documentation
29 :var __version__: version
30 :var __date__: date of the version
31 """
32
33 __author__="Lutz Gross, l.gross@uq.edu.au"
34
35
36 from types import StringType,IntType,FloatType,BooleanType,ListType,DictType
37 from sys import stdout
38 import numpy
39 import operator
40 import itertools
41 import time
42 import os
43
44 # import the 'set' module if it's not defined (python2.3/2.4 difference)
45 try:
46 set
47 except NameError:
48 from sets import Set as set
49
50 from xml.dom import minidom
51
52
53 def all(seq):
54 """
55 Returns True if no element in ``seq`` is ``None``, False otherwise.
56 """
57 for x in seq:
58 if not x:
59 return False
60 return True
61
62 def any(seq):
63 """
64 Returns True if not all elements in ``seq`` are ``None``, False otherwise.
65 """
66 for x in seq:
67 if x:
68 return True
69 return False
70
71 def importName(modulename, name):
72 """
73 Imports a named object from a module in the context of this function,
74 which means you should use fully qualified module paths.
75 Returns None on failure.
76
77 This function is from U{http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52241}
78 """
79 module = __import__(modulename, globals(), locals(), [name])
80
81 try:
82 return vars(module)[name]
83 except KeyError:
84 raise ImportError("Could not import %s from %s" % (name, modulename))
85
86 class ESySXMLParser(object):
87 """
88 Parser for an ESysXML file.
89 """
90 def __init__(self,xml, debug=False):
91 self.__dom = minidom.parseString(xml)
92 self.__linkable_object_registry= {}
93 self.__link_registry= []
94 self.__esys=self.__dom.getElementsByTagName('ESys')[0]
95 self.debug=debug
96
97 def getClassPath(self, node):
98 type = node.getAttribute("type")
99 if (node.getAttribute("module")):
100 module = node.getAttribute("module")
101 return importName(module, type)
102 else:
103 return importName("__main__", type)
104
105 def setLinks(self):
106 for obj_id, link in self.__link_registry:
107 link.target = self.__linkable_object_registry[obj_id]
108
109 def parse(self):
110 """
111 Parses EsysXML and returns the list of generating ParameterSets.
112 """
113 found=[]
114 for node in self.__esys.childNodes:
115 if isinstance(node, minidom.Element):
116 if node.tagName == 'Simulation':
117 found.append(Simulation.fromDom(self, node))
118 elif node.tagName == 'Model':
119 found.append(self.getClassPath(node).fromDom(self, node))
120 elif node.tagName == 'ParameterSet':
121 found.append(self.getClassPath(node).fromDom(self, node))
122 else:
123 raise "Invalid type, %r" % node.getAttribute("type")
124 self.setLinks()
125 return found
126
127 def registerLink(self,obj_id, link):
128 self.__link_registry.append((int(obj_id),link))
129
130 def registerLinkableObject(self,obj, node):
131 id_str=node.getAttribute('id').strip()
132 if len(id_str)>0:
133 id=int(id_str)
134 if self.__linkable_object_registry.has_key(id):
135 raise ValueError("Object id %s already exists."%id)
136 else:
137 self.__linkable_object_registry[id]=obj
138
139 def getComponent(self, node):
140 """
141 Returns a single component + rank from a simulation.
142 """
143 rank = int(node.getAttribute("rank"))
144 for n in node.childNodes:
145 if isinstance(n, minidom.Element):
146 if n.tagName == 'Simulation':
147 return (rank, Simulation.fromDom(self, n))
148 elif n.tagName == 'Model':
149 return (rank, self.getClassPath(n).fromDom(self, n))
150 elif n.tagName == 'ParameterSet':
151 return (rank, self.getClassPath(n).fromDom(self, n))
152 else:
153 raise ValueError("illegal component type %s"%n.tagName)
154 raise ValueError("cannot resolve Component")
155
156 class ESySXMLCreator(object):
157 """
158 Creates an XML Dom representation.
159 """
160 def __init__(self):
161 self.__dom=minidom.Document()
162 self.__esys =self.__dom.createElement('ESys')
163 self.__dom.appendChild(self.__esys)
164 self.__linkable_object_registry={}
165 self.__number_sequence = itertools.count(100)
166
167 def getRoot(self):
168 return self.__esys
169
170 def createElement(self,name):
171 return self.__dom.createElement(name)
172
173 def createTextNode(self,name):
174 return self.__dom.createTextNode(name)
175
176 def getElementById(self,name):
177 return self.__dom.getElementById(name)
178
179 def createDataNode(self, tagName, data):
180 """
181 ``createDataNode`` s are the building blocks of the XML documents
182 constructed in this module.
183
184 :param tagName: the associated XML tag
185 :param data: the values in the tag
186 """
187 n = self.createElement(tagName)
188 n.appendChild(self.createTextNode(str(data)))
189 return n
190
191 def getLinkableObjectId(self, obj):
192 for id, o in self.__linkable_object_registry.items():
193 if o == obj: return id
194 id =self.__number_sequence.next()
195 self.__linkable_object_registry[id]=obj
196 return id
197
198 def registerLinkableObject(self, obj, node):
199 """
200 Returns a unique object id for object ``obj``.
201 """
202 id=self.getLinkableObjectId(obj)
203 node.setAttribute('id',str(id))
204 node.setIdAttribute("id")
205
206 def includeTargets(self):
207 target_written=True
208 while target_written:
209 targetsList =self.__dom.getElementsByTagName('Target')
210 target_written=False
211 for element in targetsList:
212 targetId = int(element.firstChild.nodeValue.strip())
213 if self.getElementById(str(targetId)): continue
214 targetObj = self.__linkable_object_registry[targetId]
215 targetObj.toDom(self, self.__esys)
216 target_written=True
217
218 def toprettyxml(self):
219 self.includeTargets()
220 return self.__dom.toprettyxml()
221
222 class Link:
223 """
224 A Link makes an attribute of an object callable::
225
226 o.object()
227 o.a=8
228 l=Link(o,"a")
229 assert l()==8
230 """
231
232 def __init__(self,target,attribute=None):
233 """
234 Creates a link to the object target. If attribute is given, the link is
235 established to this attribute of the target. Otherwise the attribute is
236 undefined.
237 """
238 self.target = target
239 self.attribute = None
240 self.setAttributeName(attribute)
241
242 def getTarget(self):
243 """
244 Returns the target.
245 """
246 return self.target
247
248 def getAttributeName(self):
249 """
250 Returns the name of the attribute the link is pointing to.
251 """
252 return self.attribute
253
254 def setAttributeName(self,attribute):
255 """
256 Sets a new attribute name to be collected from the target object. The
257 target object must have the attribute with name attribute.
258 """
259 if attribute and self.target:
260 if isinstance(self.target,LinkableObject):
261 if not self.target.hasAttribute(attribute):
262 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
263 else:
264 if not hasattr(self.target,attribute):
265 raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
266 self.attribute = attribute
267
268 def hasDefinedAttributeName(self):
269 """
270 Returns true if an attribute name is set.
271 """
272 return self.attribute != None
273
274 def __repr__(self):
275 """
276 Returns a string representation of the link.
277 """
278 if self.hasDefinedAttributeName():
279 return "<Link to attribute %s of %s>" % (self.attribute, self.target)
280 else:
281 return "<Link to target %s>" % self.target
282
283 def __call__(self,name=None):
284 """
285 Returns the value of the attribute of the target object. If the
286 attribute is callable then the return value of the call is returned.
287 """
288 if name:
289 out=getattr(self.target, name)
290 else:
291 out=getattr(self.target, self.attribute)
292
293 if callable(out):
294 return out()
295 else:
296 return out
297
298 def toDom(self, esysxml, node):
299 """
300 ``toDom`` method of Link. Creates a Link node and appends it to the
301 current XML esysxml.
302 """
303 link = esysxml.createElement('Link')
304 assert (self.target != None), ("Target was none, name was %r" % self.attribute)
305 link.appendChild(esysxml.createDataNode('Target', esysxml.getLinkableObjectId(self.target)))
306 # this use of id will not work for purposes of being able to retrieve the intended
307 # target from the xml later. I need a better unique identifier.
308 assert self.attribute, "You can't xmlify a Link without a target attribute"
309 link.appendChild(esysxml.createDataNode('Attribute', self.attribute))
310 node.appendChild(link)
311
312 def fromDom(cls, esysxml, node):
313 targetid = int(node.getElementsByTagName("Target")[0].firstChild.nodeValue.strip())
314 attribute =str(node.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip())
315 l = cls(None, attribute)
316 esysxml.registerLink(targetid, l)
317 return l
318
319 fromDom = classmethod(fromDom)
320
321 class LinkableObject(object):
322 """
323 An object that allows to link its attributes to attributes of other objects
324 via a Link object. For instance::
325
326 p = LinkableObject()
327 p.x = Link(o,"name")
328 print p.x
329
330 links attribute ``x`` of ``p`` to the attribute name of object ``o``.
331
332 ``p.x`` will contain the current value of attribute ``name`` of object ``o``.
333
334 If the value of ``getattr(o, "name")`` is callable, ``p.x`` will return
335 the return value of the call.
336 """
337
338 def __init__(self, id = None, debug=False):
339 """
340 Initializes LinkableObject so that we can operate on Links.
341 """
342 self.debug = debug
343 self.__linked_attributes={}
344
345 def trace(self, msg):
346 """
347 If debugging is on, prints the message, otherwise does nothing.
348 """
349 if self.debug:
350 print "%s: %s"%(str(self),msg)
351
352 def __getattr__(self,name):
353 """
354 Returns the value of attribute name. If the value is a Link object the
355 object is called and the return value is returned.
356 """
357 out = self.getAttributeObject(name)
358 if isinstance(out,Link):
359 return out()
360 else:
361 return out
362
363 def getAttributeObject(self,name):
364 """
365 Returns the object stored for attribute ``name``.
366 """
367
368 if self.__dict__.has_key(name):
369 return self.__dict__[name]
370
371 if self.__linked_attributes.has_key(name):
372 return self.__linked_attributes[name]
373
374 if self.__class__.__dict__.has_key(name):
375 return self.__class.__dict__[name]
376
377 raise AttributeError,"No attribute %s."%name
378
379 def hasAttribute(self,name):
380 """
381 Returns True if self has attribute ``name``.
382 """
383 return self.__dict__.has_key(name) or self.__linked_attributes.has_key(name) or self.__class__.__dict__.has_key(name)
384
385 def __setattr__(self,name,value):
386 """
387 Sets the value for attribute name. If value is a Link the target
388 attribute is set to name if no attribute has been specified.
389 """
390
391 if self.__dict__.has_key(name):
392 del self.__dict__[name]
393
394 if isinstance(value,Link):
395 if not value.hasDefinedAttributeName():
396 value.setAttributeName(name)
397 self.__linked_attributes[name] = value
398
399 self.trace("attribute %s is now linked by %s."%(name,value))
400 else:
401 self.__dict__[name] = value
402
403 def __delattr__(self,name):
404 """
405 Removes the attribute ``name``.
406 """
407
408 if self.__linked_attributes.has_key[name]:
409 del self.__linked_attributes[name]
410 elif self.__dict__.has_key(name):
411 del self.__dict__[name]
412 else:
413 raise AttributeError,"No attribute %s."%name
414
415 class _ParameterIterator:
416 def __init__(self,parameterset):
417
418 self.__set=parameterset
419 self.__iter=iter(parameterset.parameters)
420
421 def next(self):
422 o=self.__iter.next()
423 return (o,self.__set.getAttributeObject(o))
424
425 def __iter__(self):
426 return self
427
428 class ParameterSet(LinkableObject):
429 """
430 A class which allows to emphasize attributes to be written and read to XML.
431
432 Leaves of an ESySParameters object can be:
433
434 - a real number
435 - an integer number
436 - a string
437 - a boolean value
438 - a ParameterSet object
439 - a Simulation object
440 - a Model object
441 - a numpy object
442 - a list of booleans
443 - any other object (not considered by writeESySXML and writeXML)
444
445 Example for how to create an ESySParameters object::
446
447 p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
448 p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
449 parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))
450
451 This can be accessed as::
452
453 parm.parm1.gamma=0.
454 parm.parm1.dim=2
455 parm.parm1.tol_v=0.001
456 parm.parm1.output_file="/tmp/u.%3.3d.dx"
457 parm.parm1.runFlag=True
458 parm.parm1.parm11.gamma1=1.
459 parm.parm1.parm11.gamma2=2.
460 parm.parm1.parm11.gamma3=3.
461 parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
462 """
463 def __init__(self, parameters=[], **kwargs):
464 """
465 Creates a ParameterSet with given parameters.
466 """
467 LinkableObject.__init__(self, **kwargs)
468 self.parameters = set()
469 self.declareParameters(parameters)
470
471 def __repr__(self):
472 return "<%s %d>"%(self.__class__.__name__,id(self))
473
474 def declareParameter(self,**parameters):
475 """
476 Declares one or more new parameters and their initial value.
477 """
478 self.declareParameters(parameters)
479
480 def declareParameters(self,parameters):
481 """
482 Declares a set of parameters. parameters can be a list, a dictionary
483 or a ParameterSet.
484 """
485 if isinstance(parameters,ListType):
486 parameters = zip(parameters, itertools.repeat(None))
487 if isinstance(parameters,DictType):
488 parameters = parameters.iteritems()
489
490 for prm, value in parameters:
491 setattr(self,prm,value)
492 self.parameters.add(prm)
493
494 def releaseParameters(self,name):
495 """
496 Removes parameter name from the parameters.
497 """
498 if self.isParameter(name):
499 self.parameters.remove(name)
500 self.trace("parameter %s has been removed."%name)
501
502 def checkLinkTargets(self, models, hash):
503 """
504 Returns a set of tuples ("<self>(<name>)", <target model>) if the
505 parameter <name> is linked to model <target model> but <target model>
506 is not in the list of models. If a parameter is linked to another
507 parameter set which is not in the hash list the parameter set is
508 checked for its models. hash gives the call history.
509 """
510 out=set()
511 for name, value in self:
512 if isinstance(value, Link):
513 m=value.getTarget()
514 if isinstance(m, Model):
515 if not m in models: out.add( (str(self)+"("+name+")",m) )
516 elif isinstance(m, ParameterSet) and not m in hash:
517 out|=set( [ (str(self)+"("+name+")."+f[0],f[1]) for f in m.checkLinkTargets(models, hash+[ self ] ) ] )
518 return out
519
520 def __iter__(self):
521 """
522 Creates an iterator over the parameter and their values.
523 """
524 return _ParameterIterator(self)
525
526 def showParameters(self):
527 """
528 Returns a description of the parameters.
529 """
530 out="{"
531 notfirst=False
532 for i,v in self:
533 if notfirst: out=out+","
534 notfirst=True
535 if isinstance(v,ParameterSet):
536 out="%s\"%s\" : %s"%(out,i,v.showParameters())
537 else:
538 out="%s\"%s\" : %s"%(out,i,v)
539 return out+"}"
540
541 def __delattr__(self,name):
542 """
543 Removes the attribute ``name``.
544 """
545 LinkableObject.__delattr__(self,name)
546 try:
547 self.releaseParameter(name)
548 except:
549 pass
550
551 def toDom(self, esysxml, node):
552 """
553 ``toDom`` method of Model class.
554 """
555 pset = esysxml.createElement('ParameterSet')
556 pset.setAttribute('type', self.__class__.__name__)
557 pset.setAttribute('module', self.__class__.__module__)
558 esysxml.registerLinkableObject(self, pset)
559 self._parametersToDom(esysxml, pset)
560 node.appendChild(pset)
561
562 def _parametersToDom(self, esysxml, node):
563 for name,value in self:
564 # convert list to numpy when possible:
565 if isinstance (value, list):
566 elem_type=-1
567 for i in value:
568 if isinstance(i, bool):
569 elem_type = max(elem_type,0)
570 elif isinstance(i, int):
571 elem_type = max(elem_type,1)
572 elif isinstance(i, float):
573 elem_type = max(elem_type,2)
574 if elem_type == 0: value = numpy.array(value,numpy.bool_)
575 if elem_type == 1: value = numpy.array(value,numpy.int_)
576 if elem_type == 2: value = numpy.array(value,numpy.float_)
577
578 param = esysxml.createElement('Parameter')
579 param.setAttribute('type', value.__class__.__name__)
580
581 param.appendChild(esysxml.createDataNode('Name', name))
582
583 val = esysxml.createElement('Value')
584 if isinstance(value,(ParameterSet,Link,DataSource)):
585 value.toDom(esysxml, val)
586 param.appendChild(val)
587 elif isinstance(value, numpy.ndarray):
588 shape = value.shape
589 if isinstance(shape, tuple):
590 size = reduce(operator.mul, shape)
591 shape = ' '.join(map(str, shape))
592 else:
593 size = shape
594 shape = str(shape)
595
596 arraytype = value.dtype.kind
597 if arraytype=='b':
598 arraytype_str="bool_"
599 elif arraytype=='i':
600 arraytype_str="int_"
601 elif arraytype=='f':
602 arraytype_str="float_"
603 elif arraytype=='c':
604 arraytype_str="complex_"
605 else:
606 arraytype_str=str(arraytype)
607 ndarrayElement = esysxml.createElement('ndarray')
608 ndarrayElement.appendChild(esysxml.createDataNode('ArrayType', arraytype_str))
609 ndarrayElement.appendChild(esysxml.createDataNode('Shape', shape))
610 ndarrayElement.appendChild(esysxml.createDataNode('Data', ' '.join(
611 [str(x) for x in numpy.reshape(value, size)])))
612 val.appendChild(ndarrayElement)
613 param.appendChild(val)
614 elif isinstance(value, list):
615 param.appendChild(esysxml.createDataNode('Value', ' '.join([str(x) for x in value]) ))
616 elif isinstance(value, (str, bool, int, float, type(None))):
617 param.appendChild(esysxml.createDataNode('Value', str(value)))
618 elif isinstance(value, dict):
619 dic = esysxml.createElement('dictionary')
620 if len(value.keys())>0:
621 dic.setAttribute('key_type', value.keys()[0].__class__.__name__)
622 dic.setAttribute('value_type', value[value.keys()[0]].__class__.__name__)
623 for k,v in value.items():
624 i=esysxml.createElement('item')
625 i.appendChild(esysxml.createDataNode('key', k))
626 i.appendChild(esysxml.createDataNode('value', v))
627 dic.appendChild(i)
628 param.appendChild(dic)
629 else:
630 raise ValueError("cannot serialize %s type to XML."%str(value.__class__))
631
632 node.appendChild(param)
633
634 def fromDom(cls, esysxml, node):
635 # Define a host of helper functions to assist us.
636 def _children(node):
637 """
638 Remove the empty nodes from the children of this node.
639 """
640 ret = []
641 for x in node.childNodes:
642 if isinstance(x, minidom.Text):
643 if x.nodeValue.strip():
644 ret.append(x)
645 else:
646 ret.append(x)
647 return ret
648
649 def _floatfromValue(esysxml, node):
650 return float(node.nodeValue.strip())
651
652 def _stringfromValue(esysxml, node):
653 return str(node.nodeValue.strip())
654
655 def _intfromValue(esysxml, node):
656 return int(node.nodeValue.strip())
657
658 def _boolfromValue(esysxml, node):
659 return _boolfromstring(node.nodeValue.strip())
660
661 def _nonefromValue(esysxml, node):
662 return None
663
664 def _ndarrayfromValue(esysxml, node):
665 for node in _children(node):
666 if node.tagName == 'ArrayType':
667 arraytype = node.firstChild.nodeValue.strip()
668 if node.tagName == 'Shape':
669 shape = node.firstChild.nodeValue.strip()
670 shape = [int(x) for x in shape.split()]
671 if node.tagName == 'Data':
672 data = node.firstChild.nodeValue.strip()
673 ndtype=getattr(numpy,arraytype)
674 if ndtype==numpy.bool_:
675 data=[(x=="True") for x in data.split()]
676 else:
677 data=[ndtype(x) for x in data.split()]
678 return numpy.reshape(numpy.array(data, dtype=ndtype),
679 shape)
680
681 def _listfromValue(esysxml, node):
682 return [x for x in node.nodeValue.split()]
683
684 def _boolfromstring(s):
685 if s == 'True':
686 return True
687 else:
688 return False
689 # Mapping from text types in the xml to methods used to process trees of that type
690 ptypemap = {"Simulation": Simulation.fromDom,
691 "Model":Model.fromDom,
692 "ParameterSet":ParameterSet.fromDom,
693 "Link":Link.fromDom,
694 "DataSource":DataSource.fromDom,
695 "float":_floatfromValue,
696 "int":_intfromValue,
697 "str":_stringfromValue,
698 "bool":_boolfromValue,
699 "list":_listfromValue,
700 "ndarray" : _ndarrayfromValue,
701 "NoneType":_nonefromValue,
702 }
703
704 parameters = {}
705 for n in _children(node):
706 ptype = n.getAttribute("type")
707 if not ptypemap.has_key(ptype):
708 raise KeyError("cannot handle parameter type %s."%ptype)
709
710 pname = pvalue = None
711 for childnode in _children(n):
712 if childnode.tagName == "Name":
713 pname = childnode.firstChild.nodeValue.strip()
714
715 if childnode.tagName == "Value":
716 nodes = _children(childnode)
717 pvalue = ptypemap[ptype](esysxml, nodes[0])
718
719 parameters[pname] = pvalue
720
721 # Create the instance of ParameterSet
722 try:
723 o = cls(debug=esysxml.debug)
724 except TypeError, inst:
725 print inst.args[0]
726 if inst.args[0]=="__init__() got an unexpected keyword argument 'debug'":
727 raise TypeError("The Model class %s __init__ needs to have argument 'debug'.")
728 else:
729 raise inst
730 o.declareParameters(parameters)
731 esysxml.registerLinkableObject(o, node)
732 return o
733
734 fromDom = classmethod(fromDom)
735
736 def writeXML(self,ostream=stdout):
737 """
738 Writes the object as an XML object into an output stream.
739 """
740 esysxml=ESySXMLCreator()
741 self.toDom(esysxml, esysxml.getRoot())
742 ostream.write(esysxml.toprettyxml())
743
744 class Model(ParameterSet):
745 """
746 A Model object represents a process marching over time until a
747 finalizing condition is fulfilled. At each time step an iterative
748 process can be performed and the time step size can be controlled. A
749 Model has the following work flow::
750
751 doInitialization()
752 while not terminateInitialIteration(): doInitializationiStep()
753 doInitialPostprocessing()
754 while not finalize():
755 dt=getSafeTimeStepSize(dt)
756 doStepPreprocessing(dt)
757 while not terminateIteration(): doStep(dt)
758 doStepPostprocessing(dt)
759 doFinalization()
760
761 where ``doInitialization``, ``finalize``, ``getSafeTimeStepSize``,
762 ``doStepPreprocessing``, ``terminateIteration``, ``doStepPostprocessing``,
763 ``doFinalization`` are methods of the particular instance of a Model. The
764 default implementations of these methods have to be overwritten by the
765 subclass implementing a Model.
766 """
767
768 UNDEF_DT=1.e300
769
770 def __init__(self,parameters=[],**kwargs):
771 """
772 Creates a model.
773
774 Just calls the parent constructor.
775 """
776 ParameterSet.__init__(self, parameters=parameters,**kwargs)
777
778 def __str__(self):
779 return "<%s %d>"%(self.__class__.__name__,id(self))
780
781
782 def setUp(self):
783 """
784 Sets up the model.
785
786 This function may be overwritten.
787 """
788 pass
789
790 def doInitialization(self):
791 """
792 Initializes the time stepping scheme. This method is not called in
793 case of a restart.
794
795 This function may be overwritten.
796 """
797 pass
798
799 def doInitialStep(self):
800 """
801 Performs an iteration step in the initialization phase. This method
802 is not called in case of a restart.
803
804 This function may be overwritten.
805 """
806 pass
807
808 def terminateInitialIteration(self):
809 """
810 Returns True if iteration at the inital phase is terminated.
811 """
812 return True
813
814 def doInitialPostprocessing(self):
815 """
816 Finalises the initialization iteration process. This method is not
817 called in case of a restart.
818
819 This function may be overwritten.
820 """
821 pass
822
823 def getSafeTimeStepSize(self,dt):
824 """
825 Returns a time step size which can be safely used.
826
827 ``dt`` gives the previously used step size.
828
829 This function may be overwritten.
830 """
831 return self.UNDEF_DT
832
833 def finalize(self):
834 """
835 Returns False if the time stepping is finalized.
836
837 This function may be overwritten.
838 """
839 return False
840
841 def doFinalization(self):
842 """
843 Finalizes the time stepping.
844
845 This function may be overwritten.
846 """
847 pass
848
849 def doStepPreprocessing(self,dt):
850 """
851 Sets up a time step of step size dt.
852
853 This function may be overwritten.
854 """
855 pass
856
857 def doStep(self,dt):
858 """
859 Executes an iteration step at a time step.
860
861 ``dt`` is the currently used time step size.
862
863 This function may be overwritten.
864 """
865 pass
866
867 def terminateIteration(self):
868 """
869 Returns True if iteration on a time step is terminated.
870 """
871 return True
872
873 def doStepPostprocessing(self,dt):
874 """
875 Finalises the time step.
876
877 dt is the currently used time step size.
878
879 This function may be overwritten.
880 """
881 pass
882
883 def toDom(self, esysxml, node):
884 """
885 ``toDom`` method of Model class.
886 """
887 pset = esysxml.createElement('Model')
888 pset.setAttribute('type', self.__class__.__name__)
889 pset.setAttribute('module', self.__class__.__module__)
890 esysxml.registerLinkableObject(self, pset)
891 node.appendChild(pset)
892 self._parametersToDom(esysxml, pset)
893
894 class Simulation(Model):
895 """
896 A Simulation object is a special Model which runs a sequence of Models.
897
898 The methods ``doInitialization``, ``finalize``, ``getSafeTimeStepSize``,
899 ``doStepPreprocessing``, ``terminateIteration``, ``doStepPostprocessing``,
900 ``doFinalization`` execute the corresponding methods of the models in
901 the simulation.
902 """
903
904 FAILED_TIME_STEPS_MAX=20
905 MAX_ITER_STEPS=50
906 MAX_CHANGE_OF_DT=2.
907
908 def __init__(self, models=[], **kwargs):
909 """
910 Initiates a simulation from a list of models.
911 """
912 super(Simulation, self).__init__(**kwargs)
913 self.declareParameter(time=0.,
914 time_step=0,
915 dt = self.UNDEF_DT)
916 for m in models:
917 if not isinstance(m, Model):
918 raise TypeError("%s is not a subclass of Model."%m)
919 self.__models=[]
920 for i in range(len(models)):
921 self[i] = models[i]
922
923
924 def __repr__(self):
925 """
926 Returns a string representation of the Simulation.
927 """
928 return "<Simulation %r>" % self.__models
929
930 def __str__(self):
931 """
932 Returns Simulation as a string.
933 """
934 return "<Simulation %d>" % id(self)
935
936 def iterModels(self):
937 """
938 Returns an iterator over the models.
939 """
940 return self.__models
941
942 def __getitem__(self,i):
943 """
944 Returns the i-th model.
945 """
946 return self.__models[i]
947
948 def __setitem__(self,i,value):
949 """
950 Sets the i-th model.
951 """
952 if not isinstance(value,Model):
953 raise ValueError,"assigned value is not a Model but instance of %s"%(value.__class__.__name__,)
954 for j in range(max(i-len(self.__models)+1,0)):
955 self.__models.append(None)
956 self.__models[i]=value
957
958 def __len__(self):
959 """
960 Returns the number of models.
961 """
962 return len(self.__models)
963
964 def getAllModels(self):
965 """
966 Returns a list of all models used in the Simulation including
967 subsimulations.
968 """
969 out=[]
970 for m in self.iterModels():
971 if isinstance(m, Simulation):
972 out+=m.getAllModels()
973 else:
974 out.append(m)
975 return list(set(out))
976
977 def checkModels(self, models, hash):
978 """
979 Returns a list of (model, parameter, target model) if the parameter
980 of model is linking to the target_model which is not in the list of
981 models.
982 """
983 out=self.checkLinkTargets(models, hash + [self])
984 for m in self.iterModels():
985 if isinstance(m, Simulation):
986 out|=m.checkModels(models, hash)
987 else:
988 out|=m.checkLinkTargets(models, hash + [self])
989 return set( [ (str(self)+"."+f[0],f[1]) for f in out ] )
990
991
992 def getSafeTimeStepSize(self,dt):
993 """
994 Returns a time step size which can be safely used by all models.
995
996 This is the minimum over the time step sizes of all models.
997 """
998 out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
999 return out
1000
1001 def setUp(self):
1002 """
1003 Performs the setup for all models.
1004 """
1005 for o in self.iterModels():
1006 o.setUp()
1007
1008 def doInitialization(self):
1009 """
1010 Initializes all models.
1011 """
1012 for o in self.iterModels():
1013 o.doInitialization()
1014
1015 def doInitialStep(self):
1016 """
1017 Performs an iteration step in the initialization step for all models.
1018 """
1019 iter=0
1020 while not self.terminateInitialIteration():
1021 if iter==0: self.trace("iteration for initialization starts")
1022 iter+=1
1023 self.trace("iteration step %d"%(iter))
1024 for o in self.iterModels():
1025 o.doInitialStep()
1026 if iter>self.MAX_ITER_STEPS:
1027 raise IterationDivergenceError("initial iteration did not converge after %s steps."%iter)
1028 self.trace("Initialization finalized after %s iteration steps."%iter)
1029
1030 def doInitialPostprocessing(self):
1031 """
1032 Finalises the initialization iteration process for all models.
1033 """
1034 for o in self.iterModels():
1035 o.doInitialPostprocessing()
1036
1037 def finalize(self):
1038 """
1039 Returns True if any of the models is to be finalized.
1040 """
1041 return any([o.finalize() for o in self.iterModels()])
1042
1043 def doFinalization(self):
1044 """
1045 Finalises the time stepping for all models.
1046 """
1047 for i in self.iterModels(): i.doFinalization()
1048 self.trace("end of time integation.")
1049
1050 def doStepPreprocessing(self,dt):
1051 """
1052 Initializes the time step for all models.
1053 """
1054 for o in self.iterModels():
1055 o.doStepPreprocessing(dt)
1056
1057 def terminateIteration(self):
1058 """
1059 Returns True if all iterations for all models are terminated.
1060 """
1061 out=all([o.terminateIteration() for o in self.iterModels()])
1062 return out
1063
1064 def terminateInitialIteration(self):
1065 """
1066 Returns True if all initial iterations for all models are terminated.
1067 """
1068 out=all([o.terminateInitialIteration() for o in self.iterModels()])
1069 return out
1070
1071 def doStepPostprocessing(self,dt):
1072 """
1073 Finalises the iteration process for all models.
1074 """
1075 for o in self.iterModels():
1076 o.doStepPostprocessing(dt)
1077 self.time_step+=1
1078 self.time+=dt
1079 self.dt=dt
1080
1081 def doStep(self,dt):
1082 """
1083 Executes the iteration step at a time step for all models::
1084
1085 self.doStepPreprocessing(dt)
1086 while not self.terminateIteration():
1087 for all models:
1088 self.doStep(dt)
1089 self.doStepPostprocessing(dt)
1090 """
1091 self.iter=0
1092 while not self.terminateIteration():
1093 if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.time_step+1,self.time+dt))
1094 self.iter+=1
1095 self.trace("iteration step %d"%(self.iter))
1096 for o in self.iterModels():
1097 o.doStep(dt)
1098 if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.time_step+1,self.time+dt))
1099
1100 def run(self,check_pointing=None):
1101 """
1102 Runs the simulation by performing essentially::
1103
1104 self.setUp()
1105 if not restart:
1106 self.doInitialization()
1107 while not self.terminateInitialIteration(): self.doInitialStep()
1108 self.doInitialPostprocessing()
1109 while not self.finalize():
1110 dt=self.getSafeTimeStepSize()
1111 self.doStepPreprocessing(dt_new)
1112 self.doStep(dt_new)
1113 self.doStepPostprocessing(dt_new)
1114 self.doFinalization()
1115
1116 If one of the models throws a ``FailedTimeStepError`` exception a
1117 new time step size is computed through getSafeTimeStepSize() and the
1118 time step is repeated.
1119
1120 If one of the models throws a ``IterationDivergenceError``
1121 exception the time step size is halved and the time step is repeated.
1122
1123 In both cases the time integration is given up after
1124 ``Simulation.FAILED_TIME_STEPS_MAX`` attempts.
1125 """
1126 # check the completeness of the models:
1127 # first a list of all the models involved in the simulation including
1128 # subsimulations:
1129 #
1130 missing=self.checkModels(self.getAllModels(), [])
1131 if len(missing)>0:
1132 msg=""
1133 for l in missing:
1134 msg+="\n\t"+str(l[1])+" at "+l[0]
1135 raise MissingLink("link targets missing in the Simulation: %s"%msg)
1136 #==============================
1137 self.setUp()
1138 if self.time_step < 1:
1139 self.doInitialization()
1140 self.doInitialStep()
1141 self.doInitialPostprocessing()
1142 while not self.finalize():
1143 step_fail_counter=0
1144 iteration_fail_counter=0
1145 if self.time_step==0:
1146 dt_new=self.getSafeTimeStepSize(self.dt)
1147 else:
1148 dt_new=min(max(self.getSafeTimeStepSize(self.dt),self.dt/self.MAX_CHANGE_OF_DT),self.dt*self.MAX_CHANGE_OF_DT)
1149 self.trace("%d. time step %e (step size %e.)" % (self.time_step+1,self.time+dt_new,dt_new))
1150 end_of_step=False
1151 while not end_of_step:
1152 end_of_step=True
1153 if not dt_new>0:
1154 raise NonPositiveStepSizeError("non-positive step size in step %d"%(self.time_step+1))
1155 try:
1156 self.doStepPreprocessing(dt_new)
1157 self.doStep(dt_new)
1158 self.doStepPostprocessing(dt_new)
1159 except IterationDivergenceError:
1160 dt_new*=0.5
1161 end_of_step=False
1162 iteration_fail_counter+=1
1163 if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
1164 raise SimulationBreakDownError("reduction of time step to achieve convergence failed after %s steps."%self.FAILED_TIME_STEPS_MAX)
1165 self.trace("Iteration failed. Time step is repeated with new step size %s."%dt_new)
1166 except FailedTimeStepError:
1167 dt_new=self.getSafeTimeStepSize(self.dt)
1168 end_of_step=False
1169 step_fail_counter+=1
1170 self.trace("Time step is repeated with new time step size %s."%dt_new)
1171 if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
1172 raise SimulationBreakDownError("Time integration is given up after %d attempts."%step_fail_counter)
1173 if not check_pointing==None:
1174 if check_pointing.doDump():
1175 self.trace("check point is created.")
1176 self.writeXML()
1177 self.doFinalization()
1178
1179
1180 def toDom(self, esysxml, node):
1181 """
1182 ``toDom`` method of Simulation class.
1183 """
1184 simulation = esysxml.createElement('Simulation')
1185 esysxml.registerLinkableObject(self, simulation)
1186 for rank, sim in enumerate(self.iterModels()):
1187 component = esysxml.createElement('Component')
1188 component.setAttribute('rank', str(rank))
1189 sim.toDom(esysxml, component)
1190 simulation.appendChild(component)
1191 node.appendChild(simulation)
1192
1193 def fromDom(cls, esysxml, node):
1194 sims = []
1195 for n in node.childNodes:
1196 if isinstance(n, minidom.Text):
1197 continue
1198 sims.append(esysxml.getComponent(n))
1199 sims.sort(_comp)
1200 sim=cls([s[1] for s in sims], debug=esysxml.debug)
1201 esysxml.registerLinkableObject(sim, node)
1202 return sim
1203
1204 fromDom = classmethod(fromDom)
1205
1206 def _comp(a,b):
1207 if a[0]<a[1]:
1208 return 1
1209 elif a[0]>a[1]:
1210 return -1
1211 else:
1212 return 0
1213
1214 class IterationDivergenceError(Exception):
1215 """
1216 Exception which is thrown if there is no convergence of the iteration
1217 process at a time step.
1218
1219 But there is a chance that a smaller step could help to reach convergence.
1220 """
1221 pass
1222
1223 class FailedTimeStepError(Exception):
1224 """
1225 Exception which is thrown if the time step fails because of a step
1226 size that has been chosen to be too large.
1227 """
1228 pass
1229
1230 class SimulationBreakDownError(Exception):
1231 """
1232 Exception which is thrown if the simulation does not manage to progress
1233 in time.
1234 """
1235 pass
1236
1237 class NonPositiveStepSizeError(Exception):
1238 """
1239 Exception which is thrown if the step size is not positive.
1240 """
1241 pass
1242
1243 class MissingLink(Exception):
1244 """
1245 Exception which is thrown when a link is missing.
1246 """
1247 pass
1248
1249 class DataSource(object):
1250 """
1251 Class for handling data sources, including local and remote files.
1252 This class is under development.
1253 """
1254
1255 def __init__(self, uri="file.ext", fileformat="unknown"):
1256 self.uri = uri
1257 self.fileformat = fileformat
1258
1259 def toDom(self, esysxml, node):
1260 """
1261 ``toDom`` method of DataSource. Creates a DataSource node and appends
1262 it to the current XML esysxml.
1263 """
1264 ds = esysxml.createElement('DataSource')
1265 ds.appendChild(esysxml.createDataNode('URI', self.uri))
1266 ds.appendChild(esysxml.createDataNode('FileFormat', self.fileformat))
1267 node.appendChild(ds)
1268
1269 def fromDom(cls, esysxml, node):
1270 uri= str(node.getElementsByTagName("URI")[0].firstChild.nodeValue.strip())
1271 fileformat= str(node.getElementsByTagName("FileFormat")[0].firstChild.nodeValue.strip())
1272 ds = cls(uri, fileformat)
1273 return ds
1274
1275 def getLocalFileName(self):
1276 return self.uri
1277
1278 fromDom = classmethod(fromDom)
1279
1280 class RestartManager(object):
1281 """
1282 A restart manager which does two things: it decides when restart files
1283 were created (when doDump returns true) and manages directories for
1284 restart files. The method getNewDumper creates a new directory and
1285 returns its name.
1286
1287 This restart manager will decide to dump restart files every dump_step
1288 calls of doDump or if more than dump_time since the last dump has
1289 elapsed. The restart manager controls two directories for dumping
1290 restart data, namely for the current and previous dump. This way the
1291 previous dump can be used for restart in the case the current dump failed.
1292
1293 :cvar SEC: unit of seconds, for instance 5*RestartManager.SEC defines 5 seconds
1294 :cvar MIN: unit of minutes, for instance 5*RestartManager.MIN defines 5 minutes
1295 :cvar H: unit of hours, for instance 5*RestartManager.H defines 5 hours
1296 :cvar D: unit of days, for instance 5*RestartManager.D defines 5 days
1297 """
1298 SEC=1.
1299 MIN=60.
1300 H=360.
1301 D=8640.
1302 def __init__(self,dump_time=1080., dump_step=None, dumper=None):
1303 """
1304 Initializes the RestartManager.
1305
1306 :param dump_time: defines the minimum time interval in SEC between two
1307 dumps. If ``None``, time is not used as criterion.
1308 :param dump_step: defines the number of calls of doDump between two
1309 dump events. If ``None``, the call counter is not
1310 used as criterion.
1311 :param dumper: defines the directory for dumping restart files.
1312 Additionally, the directories dumper+"_bkp" and
1313 dumper+"_bkp2" are used. If the directory does not
1314 exist it is created. If dumper is not present a unique
1315 directory within the current working directory is used.
1316 """
1317 self.__dump_step=dump_time
1318 self.__dump_time=dump_step
1319 self.__counter=0
1320 self.__saveMarker()
1321 if dumper == None:
1322 self.__dumper="restart"+str(os.getpid())
1323 else:
1324 self.__dumper=dumper
1325 self.__dumper_bkp=self.__dumper+"_bkp"
1326 self.__dumper_bkp2=self.__dumper+"_bkp2"
1327 self.__current_dumper=None
1328
1329 def __saveMarker(self):
1330 self.__last_restart_time=time.time()
1331 self.__last_restart_counter=self.__counter
1332
1333 def getCurrentDumper(self):
1334 """
1335 Returns the name of the currently used dumper.
1336 """
1337 return self.__current_dumper
1338
1339 def doDump(self):
1340 """
1341 Returns true if the restart should be dumped. Use ``getNewDumper`` to
1342 retrieve the directory name to be used for dumping.
1343 """
1344 if self.__dump_step == None:
1345 if self.__dump_step == None:
1346 out = False
1347 else:
1348 out = (self.__dump_step + self.__last_restart_counter) <= self.__counter
1349 else:
1350 if dump_step == None:
1351 out = (self.__last_restart_time + self.__dump_time) <= time.time()
1352 else:
1353 out = ( (self.__dump_step + self.__last_restart_counter) <= self.__counter) \
1354 or ( (self.__last_restart_time + self.__dump_time) <= time.time() )
1355 if out: self.__saveMarker()
1356 self__counter+=1
1357
1358 def getNewDumper(self):
1359 """
1360 Creates a new directory to be used for dumping and returns its name.
1361 """
1362 if os.access(self.__dumper_bkp,os.F_OK):
1363 if os.access(self.__dumper_bkp2, os.F_OK):
1364 raise RunTimeError("please remove %s."%self.__dumper_bkp2)
1365 try:
1366 os.rename(self.__dumper_bkp, self.__dumper_bkp2)
1367 except:
1368 self.__current_dumper=self.__dumper
1369 raise RunTimeError("renaming backup directory %s failed. Use %s for restart."%(self.__dumper_bkp,self.__dumper))
1370 if os.access(self.__dumper,os.F_OK):
1371 if os.access(self.__dumper_bkp, os.F_OK):
1372 raise RunTimeError("please remove %s."%self.__dumper_bkp)
1373 try:
1374 os.rename(self.__dumper, self.__dumper_bkp)
1375 except:
1376 self.__current_dumper=self.__dumper_bkp2
1377 raise RunTimeError("moving directory %s to backup failed. Use %s for restart."%(self.__dumper,self.__dumper_bkp2))
1378 try:
1379 os.mkdir(self.__dumper)
1380 except:
1381 self.__current_dumper=self.__dumper_bkp
1382 raise RunTimeError("creating a new restart directory %s failed. Use %s for restart."%(self.__dumper,self.__dumper_bkp))
1383 if os.access(self.__dumper_bkp2, os.F_OK): os.rmdir(self.__dumper_bkp2)
1384 return self.getCurrentDumper()
1385
1386
1387 # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26