/[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 1020 - (show annotations)
Mon Mar 12 10:12:36 2007 UTC (12 years, 7 months ago) by phornby
File MIME type: text/x-python
File size: 45942 byte(s)
Added explicit destructors to all Exception classes.

Fixed an ifdef in TestCase.cpp

Made the conditional definition of M_PI in LocalOps.h
depend only on M_PI being undefined.

Replace dynamically dimensioned arrays in DataFactory & DataTagged with malloc.

sort() method of list does not take a named argument
(despite the manual's claims to the contary).


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