escript/py_src/modelframe.py


########################################################
#
# Copyright (c) 2003-2010 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2010 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

"""
Environment for implementing models in escript

:var __author__: name of author
:var __copyright__: copyrights
:var __license__: licence agreement
:var __url__: url entry point on documentation
:var __version__: version
:var __date__: date of the version
"""

__author__="Lutz Gross, l.gross@uq.edu.au"


from types import StringType,IntType,FloatType,BooleanType,ListType,DictType
from sys import stdout
import numpy
import operator
import itertools
import time
import os

# import the 'set' module if it's not defined (python2.3/2.4 difference)
try:
    set
except NameError:
    from sets import Set as set

from xml.dom import minidom


def all(seq):
    """
    Returns True if no element in ``seq`` is ``None``, False otherwise.
    """
    for x in seq:
        if not x:
            return False
    return True

def any(seq):
    """
    Returns True if not all elements in ``seq`` are ``None``, False otherwise.
    """
    for x in seq:
        if x:
            return True
    return False

def importName(modulename, name):
    """
    Imports a named object from a module in the context of this function,
    which means you should use fully qualified module paths.
    Returns None on failure.

    This function is from U{http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52241}
    """
    module = __import__(modulename, globals(), locals(), [name])

    try:
        return vars(module)[name]
    except KeyError:
        raise ImportError("Could not import %s from %s" % (name, modulename))

class ESySXMLParser(object):
    """
    Parser for an ESysXML file.
    """
    def __init__(self,xml, debug=False):
        self.__dom = minidom.parseString(xml)
        self.__linkable_object_registry= {}
        self.__link_registry=  []
        self.__esys=self.__dom.getElementsByTagName('ESys')[0]
        self.debug=debug

    def getClassPath(self, node):
        type = node.getAttribute("type")
        if (node.getAttribute("module")):
            module = node.getAttribute("module")
            return importName(module, type)
        else:
            return importName("__main__", type)

    def setLinks(self):
        for obj_id, link in self.__link_registry:
            link.target = self.__linkable_object_registry[obj_id]

    def parse(self):
       """
       Parses EsysXML and returns the list of generating ParameterSets.
       """
       found=[]
       for node in self.__esys.childNodes:
           if isinstance(node, minidom.Element):
               if node.tagName == 'Simulation':
                   found.append(Simulation.fromDom(self, node))
               elif node.tagName == 'Model':
                   found.append(self.getClassPath(node).fromDom(self, node))
               elif node.tagName == 'ParameterSet':
                   found.append(self.getClassPath(node).fromDom(self, node))
               else:
                   raise "Invalid type, %r" % node.getAttribute("type")
       self.setLinks()
       return found

    def registerLink(self,obj_id, link):
        self.__link_registry.append((int(obj_id),link))

    def registerLinkableObject(self,obj, node):
        id_str=node.getAttribute('id').strip()
        if len(id_str)>0:
           id=int(id_str)
           if self.__linkable_object_registry.has_key(id):
               raise ValueError("Object id %s already exists."%id)
           else:
               self.__linkable_object_registry[id]=obj

    def getComponent(self, node):
       """
       Returns a single component + rank from a simulation.
       """
       rank  = int(node.getAttribute("rank"))
       for n in node.childNodes:
           if isinstance(n, minidom.Element):
               if n.tagName == 'Simulation':
                        return (rank, Simulation.fromDom(self, n))
               elif n.tagName == 'Model':
                        return (rank, self.getClassPath(n).fromDom(self, n))
               elif n.tagName == 'ParameterSet':
                        return (rank, self.getClassPath(n).fromDom(self, n))
               else:
                 raise ValueError("illegal component type %s"%n.tagName)
       raise ValueError("cannot resolve Component")

class ESySXMLCreator(object):
    """
    Creates an XML Dom representation.
    """
    def __init__(self):
        self.__dom=minidom.Document()
        self.__esys =self.__dom.createElement('ESys')
        self.__dom.appendChild(self.__esys)
        self.__linkable_object_registry={}
        self.__number_sequence = itertools.count(100)

    def getRoot(self):
        return self.__esys

    def createElement(self,name):
        return self.__dom.createElement(name)

    def createTextNode(self,name):
        return self.__dom.createTextNode(name)

    def getElementById(self,name):
        return self.__dom.getElementById(name)

    def createDataNode(self, tagName, data):
        """
        ``createDataNode`` s are the building blocks of the XML documents
        constructed in this module.

        :param tagName: the associated XML tag
        :param data: the values in the tag
        """
        n = self.createElement(tagName)
        n.appendChild(self.createTextNode(str(data)))
        return n

    def getLinkableObjectId(self, obj):
        for id, o in self.__linkable_object_registry.items():
            if o == obj: return id
        id =self.__number_sequence.next()
        self.__linkable_object_registry[id]=obj
        return id

    def registerLinkableObject(self, obj, node):
        """
        Returns a unique object id for object ``obj``.
        """
        id=self.getLinkableObjectId(obj)
        node.setAttribute('id',str(id))
        node.setIdAttribute("id")

    def includeTargets(self):
        target_written=True
        while target_written:
            targetsList =self.__dom.getElementsByTagName('Target')
            target_written=False
            for element in targetsList:
               targetId = int(element.firstChild.nodeValue.strip())
               if self.getElementById(str(targetId)): continue
               targetObj = self.__linkable_object_registry[targetId]
               targetObj.toDom(self, self.__esys)
               target_written=True

    def toprettyxml(self):
        self.includeTargets()
        return self.__dom.toprettyxml()

class Link:
    """
    A Link makes an attribute of an object callable::

        o.object()
        o.a=8
        l=Link(o,"a")
        assert l()==8
    """

    def __init__(self,target,attribute=None):
        """
        Creates a link to the object target. If attribute is given, the link is
        established to this attribute of the target. Otherwise the attribute is
        undefined.
        """
        self.target = target
        self.attribute = None
        self.setAttributeName(attribute)

    def getTarget(self):
        """
        Returns the target.
        """
        return self.target

    def getAttributeName(self):
        """
        Returns the name of the attribute the link is pointing to.
        """
        return self.attribute

    def setAttributeName(self,attribute):
        """
        Sets a new attribute name to be collected from the target object. The
        target object must have the attribute with name attribute.
        """
        if attribute and self.target:
            if isinstance(self.target,LinkableObject):
               if not self.target.hasAttribute(attribute):
                  raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
            else:
               if not hasattr(self.target,attribute):
                  raise AttributeError("%s: target %s has no attribute %s."%(self, self.target, attribute))
        self.attribute = attribute

    def hasDefinedAttributeName(self):
        """
        Returns true if an attribute name is set.
        """
        return self.attribute != None

    def __repr__(self):
        """
        Returns a string representation of the link.
        """
        if self.hasDefinedAttributeName():
            return "<Link to attribute %s of %s>" % (self.attribute, self.target)
        else:
            return "<Link to target %s>" % self.target

    def __call__(self,name=None):
        """
        Returns the value of the attribute of the target object. If the
        attribute is callable then the return value of the call is returned.
        """
        if name:
            out=getattr(self.target, name)
        else:
            out=getattr(self.target, self.attribute)

        if callable(out):
            return out()
        else:
            return out

    def toDom(self, esysxml, node):
        """
        ``toDom`` method of Link. Creates a Link node and appends it to the
        current XML esysxml.
        """
        link = esysxml.createElement('Link')
        assert (self.target != None), ("Target was none, name was %r" % self.attribute)
        link.appendChild(esysxml.createDataNode('Target', esysxml.getLinkableObjectId(self.target)))
        # this use of id will not work for purposes of being able to retrieve the intended
        # target from the xml later. I need a better unique identifier.
        assert self.attribute, "You can't xmlify a Link without a target attribute"
        link.appendChild(esysxml.createDataNode('Attribute', self.attribute))
        node.appendChild(link)

    def fromDom(cls, esysxml, node):
        targetid = int(node.getElementsByTagName("Target")[0].firstChild.nodeValue.strip())
        attribute =str(node.getElementsByTagName("Attribute")[0].firstChild.nodeValue.strip())
        l = cls(None, attribute)
        esysxml.registerLink(targetid, l)
        return l

    fromDom = classmethod(fromDom)

class LinkableObject(object):
    """
    An object that allows to link its attributes to attributes of other objects
    via a Link object. For instance::

        p = LinkableObject()
        p.x = Link(o,"name")
        print p.x

    links attribute ``x`` of ``p`` to the attribute name of object ``o``.

    ``p.x`` will contain the current value of attribute ``name`` of object ``o``.

    If the value of ``getattr(o, "name")`` is callable, ``p.x`` will return
    the return value of the call.
    """

    def __init__(self, id = None, debug=False):
        """
        Initializes LinkableObject so that we can operate on Links.
        """
        self.debug = debug
        self.__linked_attributes={}

    def trace(self, msg):
        """
        If debugging is on, prints the message, otherwise does nothing.
        """
        if self.debug:
            print "%s: %s"%(str(self),msg)

    def __getattr__(self,name):
        """
        Returns the value of attribute name. If the value is a Link object the
        object is called and the return value is returned.
        """
        out = self.getAttributeObject(name)
        if isinstance(out,Link):
            return out()
        else:
            return out

    def getAttributeObject(self,name):
        """
        Returns the object stored for attribute ``name``.
        """

        if self.__dict__.has_key(name):
            return self.__dict__[name]

        if self.__linked_attributes.has_key(name):
            return self.__linked_attributes[name]

        if self.__class__.__dict__.has_key(name):
            return self.__class.__dict__[name]

        raise AttributeError,"No attribute %s."%name

    def hasAttribute(self,name):
        """
        Returns True if self has attribute ``name``.
        """
        return self.__dict__.has_key(name) or self.__linked_attributes.has_key(name) or  self.__class__.__dict__.has_key(name)

    def __setattr__(self,name,value):
        """
        Sets the value for attribute name. If value is a Link the target
        attribute is set to name if no attribute has been specified.
        """

        if self.__dict__.has_key(name):
            del self.__dict__[name]

        if isinstance(value,Link):
            if not value.hasDefinedAttributeName():
                value.setAttributeName(name)
            self.__linked_attributes[name] = value

            self.trace("attribute %s is now linked by %s."%(name,value))
        else:
            self.__dict__[name] = value

    def __delattr__(self,name):
        """
        Removes the attribute ``name``.
        """

        if self.__linked_attributes.has_key[name]:
            del self.__linked_attributes[name]
        elif self.__dict__.has_key(name):
            del self.__dict__[name]
        else:
            raise AttributeError,"No attribute %s."%name

class _ParameterIterator:
    def __init__(self,parameterset):

        self.__set=parameterset
        self.__iter=iter(parameterset.parameters)

    def next(self):
        o=self.__iter.next()
        return (o,self.__set.getAttributeObject(o))

    def __iter__(self):
        return self

class ParameterSet(LinkableObject):
    """
    A class which allows to emphasize attributes to be written and read to XML.

    Leaves of an ESySParameters object can be:

        - a real number
        - an integer number
        - a string
        - a boolean value
        - a ParameterSet object
        - a Simulation object
        - a Model object
        - a numpy object
        - a list of booleans
        - any other object (not considered by writeESySXML and writeXML)

    Example for how to create an ESySParameters object::

        p11=ParameterSet(gamma1=1.,gamma2=2.,gamma3=3.)
        p1=ParameterSet(dim=2,tol_v=0.001,output_file="/tmp/u.%3.3d.dx",runFlag=True,parm11=p11)
        parm=ParameterSet(parm1=p1,parm2=ParameterSet(alpha=Link(p11,"gamma1")))

    This can be accessed as::

        parm.parm1.gamma=0.
        parm.parm1.dim=2
        parm.parm1.tol_v=0.001
        parm.parm1.output_file="/tmp/u.%3.3d.dx"
        parm.parm1.runFlag=True
        parm.parm1.parm11.gamma1=1.
        parm.parm1.parm11.gamma2=2.
        parm.parm1.parm11.gamma3=3.
        parm.parm2.alpha=1. (value of parm.parm1.parm11.gamma1)
    """
    def __init__(self, parameters=[], **kwargs):
        """
        Creates a ParameterSet with given parameters.
        """
        LinkableObject.__init__(self, **kwargs)
        self.parameters = set()
        self.declareParameters(parameters)

    def __repr__(self):
        return "<%s %d>"%(self.__class__.__name__,id(self))

    def declareParameter(self,**parameters):
        """
        Declares one or more new parameters and their initial value.
        """
        self.declareParameters(parameters)

    def declareParameters(self,parameters):
        """
        Declares a set of parameters. parameters can be a list, a dictionary
        or a ParameterSet.
        """
        if isinstance(parameters,ListType):
            parameters = zip(parameters, itertools.repeat(None))
        if isinstance(parameters,DictType):
            parameters = parameters.iteritems()

        for prm, value in parameters:
            setattr(self,prm,value)
            self.parameters.add(prm)

    def releaseParameters(self,name):
        """
        Removes parameter name from the parameters.
        """
        if self.isParameter(name):
            self.parameters.remove(name)
            self.trace("parameter %s has been removed."%name)

    def checkLinkTargets(self, models, hash):
        """
        Returns a set of tuples ("<self>(<name>)", <target model>) if the
        parameter <name> is linked to model <target model> but <target model>
        is not in the list of models. If a parameter is linked to another
        parameter set which is not in the hash list the parameter set is
        checked for its models. hash gives the call history.
        """
        out=set()
        for name, value in self:
            if isinstance(value, Link):
               m=value.getTarget()
               if isinstance(m, Model):
                   if not m in models: out.add( (str(self)+"("+name+")",m) )
               elif isinstance(m, ParameterSet) and not m in hash:
                     out|=set( [ (str(self)+"("+name+")."+f[0],f[1]) for f in m.checkLinkTargets(models, hash+[ self ] ) ] )
        return out

    def __iter__(self):
        """
        Creates an iterator over the parameter and their values.
        """
        return _ParameterIterator(self)

    def showParameters(self):
        """
        Returns a description of the parameters.
        """
        out="{"
        notfirst=False
        for i,v in self:
            if notfirst: out=out+","
            notfirst=True
            if isinstance(v,ParameterSet):
                out="%s\"%s\" : %s"%(out,i,v.showParameters())
            else:
                out="%s\"%s\" : %s"%(out,i,v)
        return out+"}"

    def __delattr__(self,name):
        """
        Removes the attribute ``name``.
        """
        LinkableObject.__delattr__(self,name)
        try:
            self.releaseParameter(name)
        except:
            pass

    def toDom(self, esysxml, node):
        """
        ``toDom`` method of Model class.
        """
        pset = esysxml.createElement('ParameterSet')
        pset.setAttribute('type', self.__class__.__name__)
        pset.setAttribute('module', self.__class__.__module__)
        esysxml.registerLinkableObject(self, pset)
        self._parametersToDom(esysxml, pset)
        node.appendChild(pset)

    def _parametersToDom(self, esysxml, node):
        for name,value in self:
            # convert list to numpy when possible:
            if isinstance (value, list):
                elem_type=-1
                for i in value:
                    if isinstance(i, bool):
                        elem_type = max(elem_type,0)
                    elif isinstance(i, int):
                        elem_type = max(elem_type,1)
                    elif isinstance(i, float):
                        elem_type = max(elem_type,2)
                if elem_type == 0: value = numpy.array(value,numpy.bool_)
                if elem_type == 1: value = numpy.array(value,numpy.int_)
                if elem_type == 2: value = numpy.array(value,numpy.float_)

            param = esysxml.createElement('Parameter')
            param.setAttribute('type', value.__class__.__name__)

            param.appendChild(esysxml.createDataNode('Name', name))

            val = esysxml.createElement('Value')
            if isinstance(value,(ParameterSet,Link,DataSource)):
                value.toDom(esysxml, val)
                param.appendChild(val)
            elif isinstance(value, numpy.ndarray):
                shape = value.shape
                if isinstance(shape, tuple):
                    size = reduce(operator.mul, shape)
                    shape = ' '.join(map(str, shape))
                else:
                    size = shape
                    shape = str(shape)

                arraytype = value.dtype.kind
                if arraytype=='b':
                      arraytype_str="bool_"
                elif arraytype=='i':
                      arraytype_str="int_"
                elif arraytype=='f':
                      arraytype_str="float_"
                elif arraytype=='c':
                      arraytype_str="complex_"
                else:
                      arraytype_str=str(arraytype)
                ndarrayElement = esysxml.createElement('ndarray')
                ndarrayElement.appendChild(esysxml.createDataNode('ArrayType', arraytype_str))
                ndarrayElement.appendChild(esysxml.createDataNode('Shape', shape))
                ndarrayElement.appendChild(esysxml.createDataNode('Data', ' '.join(
                    [str(x) for x in numpy.reshape(value, size)])))
                val.appendChild(ndarrayElement)
                param.appendChild(val)
            elif isinstance(value, list):
                param.appendChild(esysxml.createDataNode('Value', ' '.join([str(x) for x in value]) ))
            elif isinstance(value, (str, bool, int, float, type(None))):
                param.appendChild(esysxml.createDataNode('Value', str(value)))
            elif isinstance(value, dict):
                 dic = esysxml.createElement('dictionary')
                 if len(value.keys())>0:
                     dic.setAttribute('key_type', value.keys()[0].__class__.__name__)
                     dic.setAttribute('value_type', value[value.keys()[0]].__class__.__name__)
                 for k,v in value.items():
                    i=esysxml.createElement('item')
                    i.appendChild(esysxml.createDataNode('key', k))
                    i.appendChild(esysxml.createDataNode('value', v))
                    dic.appendChild(i)
                 param.appendChild(dic)
            else:
                raise ValueError("cannot serialize %s type to XML."%str(value.__class__))

            node.appendChild(param)

    def fromDom(cls, esysxml, node):
        # Define a host of helper functions to assist us.
        def _children(node):
            """
            Remove the empty nodes from the children of this node.
            """
            ret = []
            for x in node.childNodes:
                if isinstance(x, minidom.Text):
                    if x.nodeValue.strip():
                        ret.append(x)
                else:
                    ret.append(x)
            return ret

        def _floatfromValue(esysxml, node):
            return float(node.nodeValue.strip())

        def _stringfromValue(esysxml, node):
            return str(node.nodeValue.strip())

        def _intfromValue(esysxml, node):
            return int(node.nodeValue.strip())

        def _boolfromValue(esysxml, node):
            return _boolfromstring(node.nodeValue.strip())

        def _nonefromValue(esysxml, node):
            return None

        def _ndarrayfromValue(esysxml, node):
            for node in _children(node):
                if node.tagName == 'ArrayType':
                    arraytype = node.firstChild.nodeValue.strip()
                if node.tagName == 'Shape':
                    shape = node.firstChild.nodeValue.strip()
                    shape = [int(x) for x in shape.split()]
                if node.tagName == 'Data':
                    data = node.firstChild.nodeValue.strip()
        ndtype=getattr(numpy,arraytype)
        if ndtype==numpy.bool_:
        data=[(x=="True") for x in data.split()]
        else:
        data=[ndtype(x) for x in data.split()]
            return numpy.reshape(numpy.array(data, dtype=ndtype),
                                    shape)

        def _listfromValue(esysxml, node):
            return [x for x in node.nodeValue.split()]

        def _boolfromstring(s):
            if s == 'True':
                return True
            else:
                return False
        # Mapping from text types in the xml to methods used to process trees of that type
        ptypemap = {"Simulation": Simulation.fromDom,
                    "Model":Model.fromDom,
                    "ParameterSet":ParameterSet.fromDom,
                    "Link":Link.fromDom,
                    "DataSource":DataSource.fromDom,
                    "float":_floatfromValue,
                    "int":_intfromValue,
                    "str":_stringfromValue,
                    "bool":_boolfromValue,
                    "list":_listfromValue,
            "ndarray" : _ndarrayfromValue,
                    "NoneType":_nonefromValue,
                    }

        parameters = {}
        for n in _children(node):
            ptype = n.getAttribute("type")
            if not ptypemap.has_key(ptype):
               raise KeyError("cannot handle parameter type %s."%ptype)

            pname = pvalue = None
            for childnode in _children(n):
                if childnode.tagName == "Name":
                    pname = childnode.firstChild.nodeValue.strip()

                if childnode.tagName == "Value":
                    nodes = _children(childnode)
                    pvalue = ptypemap[ptype](esysxml, nodes[0])

            parameters[pname] = pvalue

        # Create the instance of ParameterSet
        try:
           o = cls(debug=esysxml.debug)
        except TypeError, inst:
           print inst.args[0]
           if inst.args[0]=="__init__() got an unexpected keyword argument 'debug'":
              raise TypeError("The Model class %s __init__ needs to have argument 'debug'.")
           else:
              raise inst
        o.declareParameters(parameters)
        esysxml.registerLinkableObject(o, node)
        return o

    fromDom = classmethod(fromDom)

    def writeXML(self,ostream=stdout):
        """
        Writes the object as an XML object into an output stream.
        """
        esysxml=ESySXMLCreator()
        self.toDom(esysxml, esysxml.getRoot())
        ostream.write(esysxml.toprettyxml())

class Model(ParameterSet):
    """
    A Model object represents a process marching over time until a
    finalizing condition is fulfilled. At each time step an iterative
    process can be performed and the time step size can be controlled. A
    Model has the following work flow::

        doInitialization()
        while not terminateInitialIteration(): doInitializationiStep()
        doInitialPostprocessing()
        while not finalize():
            dt=getSafeTimeStepSize(dt)
            doStepPreprocessing(dt)
            while not terminateIteration(): doStep(dt)
            doStepPostprocessing(dt)
        doFinalization()

    where ``doInitialization``, ``finalize``, ``getSafeTimeStepSize``,
    ``doStepPreprocessing``, ``terminateIteration``, ``doStepPostprocessing``,
    ``doFinalization`` are methods of the particular instance of a Model. The
    default implementations of these methods have to be overwritten by the
    subclass implementing a Model.
    """

    UNDEF_DT=1.e300

    def __init__(self,parameters=[],**kwargs):
        """
        Creates a model.

        Just calls the parent constructor.
        """
        ParameterSet.__init__(self, parameters=parameters,**kwargs)

    def __str__(self):
       return "<%s %d>"%(self.__class__.__name__,id(self))


    def setUp(self):
        """
        Sets up the model.

        This function may be overwritten.
        """
        pass

    def doInitialization(self):
        """
        Initializes the time stepping scheme. This method is not called in
        case of a restart.

        This function may be overwritten.
        """
        pass

    def doInitialStep(self):
        """
        Performs an iteration step in the initialization phase. This method
        is not called in case of a restart.

        This function may be overwritten.
        """
        pass

    def terminateInitialIteration(self):
        """
        Returns True if iteration at the inital phase is terminated.
        """
        return True

    def doInitialPostprocessing(self):
        """
        Finalises the initialization iteration process. This method is not
        called in case of a restart.

        This function may be overwritten.
        """
        pass

    def getSafeTimeStepSize(self,dt):
        """
        Returns a time step size which can be safely used.

        ``dt`` gives the previously used step size.

        This function may be overwritten.
        """
        return self.UNDEF_DT

    def finalize(self):
        """
        Returns False if the time stepping is finalized.

        This function may be overwritten.
        """
        return False

    def doFinalization(self):
        """
        Finalizes the time stepping.

        This function may be overwritten.
        """
        pass

    def doStepPreprocessing(self,dt):
        """
        Sets up a time step of step size dt.

        This function may be overwritten.
        """
        pass

    def doStep(self,dt):
        """
        Executes an iteration step at a time step.

        ``dt`` is the currently used time step size.

        This function may be overwritten.
        """
        pass

    def terminateIteration(self):
        """
        Returns True if iteration on a time step is terminated.
        """
        return True

    def doStepPostprocessing(self,dt):
        """
        Finalises the time step.

        dt is the currently used time step size.

        This function may be overwritten.
        """
        pass

    def toDom(self, esysxml, node):
        """
        ``toDom`` method of Model class.
        """
        pset = esysxml.createElement('Model')
        pset.setAttribute('type', self.__class__.__name__)
        pset.setAttribute('module', self.__class__.__module__)
        esysxml.registerLinkableObject(self, pset)
        node.appendChild(pset)
        self._parametersToDom(esysxml, pset)

class Simulation(Model):
    """
    A Simulation object is a special Model which runs a sequence of Models.

    The methods ``doInitialization``, ``finalize``, ``getSafeTimeStepSize``,
    ``doStepPreprocessing``, ``terminateIteration``, ``doStepPostprocessing``,
    ``doFinalization`` execute the corresponding methods of the models in
    the simulation.
    """

    FAILED_TIME_STEPS_MAX=20
    MAX_ITER_STEPS=50
    MAX_CHANGE_OF_DT=2.

    def __init__(self, models=[], **kwargs):
        """
        Initiates a simulation from a list of models.
        """
        super(Simulation, self).__init__(**kwargs)
        self.declareParameter(time=0.,
                              time_step=0,
                              dt = self.UNDEF_DT)
        for m in models:
            if not isinstance(m, Model):
                 raise TypeError("%s is not a subclass of Model."%m)
        self.__models=[]
        for i in range(len(models)):
            self[i] = models[i]


    def __repr__(self):
        """
        Returns a string representation of the Simulation.
        """
        return "<Simulation %r>" % self.__models

    def __str__(self):
        """
        Returns Simulation as a string.
        """
        return "<Simulation %d>" % id(self)

    def iterModels(self):
        """
        Returns an iterator over the models.
        """
        return self.__models

    def __getitem__(self,i):
        """
        Returns the i-th model.
        """
        return self.__models[i]

    def __setitem__(self,i,value):
        """
        Sets the i-th model.
        """
        if not isinstance(value,Model):
            raise ValueError,"assigned value is not a Model but instance of %s"%(value.__class__.__name__,)
        for j in range(max(i-len(self.__models)+1,0)):
            self.__models.append(None)
        self.__models[i]=value

    def __len__(self):
        """
        Returns the number of models.
        """
        return len(self.__models)

    def getAllModels(self):
        """
        Returns a list of all models used in the Simulation including
        subsimulations.
        """
        out=[]
        for m in self.iterModels():
            if isinstance(m, Simulation):
               out+=m.getAllModels()
            else:
               out.append(m)
        return list(set(out))

    def checkModels(self, models, hash):
        """
        Returns a list of (model, parameter, target model) if the parameter
        of model is linking to the target_model which is not in the list of
        models.
        """
        out=self.checkLinkTargets(models, hash + [self])
        for m in self.iterModels():
            if isinstance(m, Simulation):
                 out|=m.checkModels(models, hash)
            else:
                 out|=m.checkLinkTargets(models, hash + [self])
        return set( [ (str(self)+"."+f[0],f[1]) for f in out ] )


    def getSafeTimeStepSize(self,dt):
        """
        Returns a time step size which can be safely used by all models.

        This is the minimum over the time step sizes of all models.
        """
        out=min([o.getSafeTimeStepSize(dt) for o in self.iterModels()])
        return out

    def setUp(self):
        """
        Performs the setup for all models.
        """
        for o in self.iterModels():
             o.setUp()

    def doInitialization(self):
        """
        Initializes all models.
        """
        for o in self.iterModels():
             o.doInitialization()

    def doInitialStep(self):
        """
        Performs an iteration step in the initialization step for all models.
        """
        iter=0
        while not self.terminateInitialIteration():
            if iter==0: self.trace("iteration for initialization starts")
            iter+=1
            self.trace("iteration step %d"%(iter))
            for o in self.iterModels():
                 o.doInitialStep()
            if iter>self.MAX_ITER_STEPS:
                 raise IterationDivergenceError("initial iteration did not converge after %s steps."%iter)
        self.trace("Initialization finalized after %s iteration steps."%iter)

    def doInitialPostprocessing(self):
        """
        Finalises the initialization iteration process for all models.
        """
        for o in self.iterModels():
            o.doInitialPostprocessing()

    def finalize(self):
        """
        Returns True if any of the models is to be finalized.
        """
        return any([o.finalize() for o in self.iterModels()])

    def doFinalization(self):
        """
        Finalises the time stepping for all models.
        """
        for i in self.iterModels(): i.doFinalization()
        self.trace("end of time integation.")

    def doStepPreprocessing(self,dt):
        """
        Initializes the time step for all models.
        """
        for o in self.iterModels():
            o.doStepPreprocessing(dt)

    def terminateIteration(self):
        """
        Returns True if all iterations for all models are terminated.
        """
        out=all([o.terminateIteration() for o in self.iterModels()])
        return out

    def terminateInitialIteration(self):
        """
        Returns True if all initial iterations for all models are terminated.
        """
        out=all([o.terminateInitialIteration() for o in self.iterModels()])
        return out

    def doStepPostprocessing(self,dt):
        """
        Finalises the iteration process for all models.
        """
        for o in self.iterModels():
            o.doStepPostprocessing(dt)
        self.time_step+=1
        self.time+=dt
        self.dt=dt

    def doStep(self,dt):
        """
        Executes the iteration step at a time step for all models::

            self.doStepPreprocessing(dt)
            while not self.terminateIteration():
                for all models:
                self.doStep(dt)
                self.doStepPostprocessing(dt)
        """
        self.iter=0
        while not self.terminateIteration():
            if self.iter==0: self.trace("iteration at %d-th time step %e starts"%(self.time_step+1,self.time+dt))
            self.iter+=1
            self.trace("iteration step %d"%(self.iter))
            for o in self.iterModels():
                  o.doStep(dt)
        if self.iter>0: self.trace("iteration at %d-th time step %e finalized."%(self.time_step+1,self.time+dt))

    def run(self,check_pointing=None):
        """
        Runs the simulation by performing essentially::

            self.setUp()
            if not restart:
                self.doInitialization()
                while not self.terminateInitialIteration(): self.doInitialStep()
                self.doInitialPostprocessing()
            while not self.finalize():
                dt=self.getSafeTimeStepSize()
                self.doStepPreprocessing(dt_new)
                self.doStep(dt_new)
                self.doStepPostprocessing(dt_new)
            self.doFinalization()

        If one of the models throws a ``FailedTimeStepError`` exception a
        new time step size is computed through getSafeTimeStepSize() and the
        time step is repeated.

        If one of the models throws a ``IterationDivergenceError``
        exception the time step size is halved and the time step is repeated.

        In both cases the time integration is given up after
        ``Simulation.FAILED_TIME_STEPS_MAX`` attempts.
        """
        # check the completeness of the models:
        # first a list of all the models involved in the simulation including
        # subsimulations:
        #
        missing=self.checkModels(self.getAllModels(), [])
        if len(missing)>0:
            msg=""
            for l in missing:
                 msg+="\n\t"+str(l[1])+" at "+l[0]
            raise MissingLink("link targets missing in the Simulation: %s"%msg)
        #==============================
        self.setUp()
        if self.time_step < 1:
           self.doInitialization()
           self.doInitialStep()
           self.doInitialPostprocessing()
        while not self.finalize():
            step_fail_counter=0
            iteration_fail_counter=0
            if self.time_step==0:
                dt_new=self.getSafeTimeStepSize(self.dt)
            else:
                dt_new=min(max(self.getSafeTimeStepSize(self.dt),self.dt/self.MAX_CHANGE_OF_DT),self.dt*self.MAX_CHANGE_OF_DT)
            self.trace("%d. time step %e (step size %e.)" % (self.time_step+1,self.time+dt_new,dt_new))
            end_of_step=False
            while not end_of_step:
               end_of_step=True
               if not dt_new>0:
                  raise NonPositiveStepSizeError("non-positive step size in step %d"%(self.time_step+1))
               try:
                  self.doStepPreprocessing(dt_new)
                  self.doStep(dt_new)
                  self.doStepPostprocessing(dt_new)
               except IterationDivergenceError:
                  dt_new*=0.5
                  end_of_step=False
                  iteration_fail_counter+=1
                  if iteration_fail_counter>self.FAILED_TIME_STEPS_MAX:
                           raise SimulationBreakDownError("reduction of time step to achieve convergence failed after %s steps."%self.FAILED_TIME_STEPS_MAX)
                  self.trace("Iteration failed. Time step is repeated with new step size %s."%dt_new)
               except FailedTimeStepError:
                  dt_new=self.getSafeTimeStepSize(self.dt)
                  end_of_step=False
                  step_fail_counter+=1
                  self.trace("Time step is repeated with new time step size %s."%dt_new)
                  if step_fail_counter>self.FAILED_TIME_STEPS_MAX:
                        raise SimulationBreakDownError("Time integration is given up after %d attempts."%step_fail_counter)
            if not check_pointing==None:
               if check_pointing.doDump():
                    self.trace("check point is created.")
                    self.writeXML()
        self.doFinalization()


    def toDom(self, esysxml, node):
        """
        ``toDom`` method of Simulation class.
        """
        simulation = esysxml.createElement('Simulation')
        esysxml.registerLinkableObject(self, simulation)
        for rank, sim in enumerate(self.iterModels()):
            component = esysxml.createElement('Component')
            component.setAttribute('rank', str(rank))
            sim.toDom(esysxml, component)
            simulation.appendChild(component)
        node.appendChild(simulation)

    def fromDom(cls, esysxml, node):
        sims = []
        for n in node.childNodes:
            if isinstance(n, minidom.Text):
                continue
            sims.append(esysxml.getComponent(n))
        sims.sort(_comp)
        sim=cls([s[1] for s in sims], debug=esysxml.debug)
        esysxml.registerLinkableObject(sim, node)
        return sim

    fromDom = classmethod(fromDom)

def _comp(a,b):
    if a[0]<a[1]:
        return 1
    elif a[0]>a[1]:
        return -1
    else:
        return 0

class IterationDivergenceError(Exception):
    """
    Exception which is thrown if there is no convergence of the iteration
    process at a time step.

    But there is a chance that a smaller step could help to reach convergence.
    """
    pass

class FailedTimeStepError(Exception):
    """
    Exception which is thrown if the time step fails because of a step
    size that has been chosen to be too large.
    """
    pass

class SimulationBreakDownError(Exception):
    """
    Exception which is thrown if the simulation does not manage to progress
    in time.
    """
    pass

class NonPositiveStepSizeError(Exception):
    """
    Exception which is thrown if the step size is not positive.
    """
    pass

class MissingLink(Exception):
    """
    Exception which is thrown when a link is missing.
    """
    pass

class DataSource(object):
    """
    Class for handling data sources, including local and remote files.
    This class is under development.
    """

    def __init__(self, uri="file.ext", fileformat="unknown"):
        self.uri = uri
        self.fileformat = fileformat

    def toDom(self, esysxml, node):
        """
        ``toDom`` method of DataSource. Creates a DataSource node and appends
        it to the current XML esysxml.
        """
        ds = esysxml.createElement('DataSource')
        ds.appendChild(esysxml.createDataNode('URI', self.uri))
        ds.appendChild(esysxml.createDataNode('FileFormat', self.fileformat))
        node.appendChild(ds)

    def fromDom(cls, esysxml, node):
        uri= str(node.getElementsByTagName("URI")[0].firstChild.nodeValue.strip())
        fileformat= str(node.getElementsByTagName("FileFormat")[0].firstChild.nodeValue.strip())
        ds = cls(uri, fileformat)
        return ds

    def getLocalFileName(self):
        return self.uri

    fromDom = classmethod(fromDom)

class RestartManager(object):
     """
     A restart manager which does two things: it decides when restart files
     were created (when doDump returns true) and manages directories for
     restart files. The method getNewDumper creates a new directory and
     returns its name.

     This restart manager will decide to dump restart files every dump_step
     calls of doDump or if more than dump_time since the last dump has
     elapsed. The restart manager controls two directories for dumping
     restart data, namely for the current and previous dump. This way the
     previous dump can be used for restart in the case the current dump failed.

     :cvar SEC: unit of seconds, for instance 5*RestartManager.SEC defines 5 seconds
     :cvar MIN: unit of minutes, for instance 5*RestartManager.MIN defines 5 minutes
     :cvar H: unit of hours, for instance 5*RestartManager.H defines 5 hours
     :cvar D: unit of days, for instance 5*RestartManager.D defines 5 days
     """
     SEC=1.
     MIN=60.
     H=360.
     D=8640.
     def __init__(self,dump_time=1080., dump_step=None, dumper=None):
         """
         Initializes the RestartManager.

         :param dump_time: defines the minimum time interval in SEC between two
                           dumps. If ``None``, time is not used as criterion.
         :param dump_step: defines the number of calls of doDump between two
                           dump events. If ``None``, the call counter is not
                           used as criterion.
         :param dumper: defines the directory for dumping restart files.
                        Additionally, the directories dumper+"_bkp" and
                        dumper+"_bkp2" are used. If the directory does not
                        exist it is created. If dumper is not present a unique
                        directory within the current working directory is used.
         """
         self.__dump_step=dump_time
         self.__dump_time=dump_step
         self.__counter=0
         self.__saveMarker()
         if dumper == None:
            self.__dumper="restart"+str(os.getpid())
         else:
            self.__dumper=dumper
         self.__dumper_bkp=self.__dumper+"_bkp"
         self.__dumper_bkp2=self.__dumper+"_bkp2"
         self.__current_dumper=None

     def __saveMarker(self):
         self.__last_restart_time=time.time()
         self.__last_restart_counter=self.__counter

     def getCurrentDumper(self):
         """
         Returns the name of the currently used dumper.
         """
         return self.__current_dumper

     def doDump(self):
        """
        Returns true if the restart should be dumped. Use ``getNewDumper`` to
        retrieve the directory name to be used for dumping.
        """
        if self.__dump_step == None:
           if self.__dump_step == None:
              out = False
           else:
              out = (self.__dump_step + self.__last_restart_counter) <= self.__counter
        else:
           if dump_step == None:
              out = (self.__last_restart_time + self.__dump_time) <= time.time()
           else:
              out = ( (self.__dump_step + self.__last_restart_counter) <= self.__counter)  \
                    or ( (self.__last_restart_time + self.__dump_time) <= time.time() )
        if out: self.__saveMarker()
        self__counter+=1

     def getNewDumper(self):
       """
       Creates a new directory to be used for dumping and returns its name.
       """
       if os.access(self.__dumper_bkp,os.F_OK):
          if os.access(self.__dumper_bkp2, os.F_OK):
             raise RunTimeError("please remove %s."%self.__dumper_bkp2)
          try:
             os.rename(self.__dumper_bkp, self.__dumper_bkp2)
          except:
             self.__current_dumper=self.__dumper
             raise RunTimeError("renaming backup directory %s failed. Use %s for restart."%(self.__dumper_bkp,self.__dumper))
       if os.access(self.__dumper,os.F_OK):
          if os.access(self.__dumper_bkp, os.F_OK):
             raise RunTimeError("please remove %s."%self.__dumper_bkp)
          try:
             os.rename(self.__dumper, self.__dumper_bkp)
          except:
             self.__current_dumper=self.__dumper_bkp2
             raise RunTimeError("moving directory %s to backup failed. Use %s for restart."%(self.__dumper,self.__dumper_bkp2))
       try:
          os.mkdir(self.__dumper)
       except:
          self.__current_dumper=self.__dumper_bkp
          raise RunTimeError("creating a new restart directory %s failed. Use %s for restart."%(self.__dumper,self.__dumper_bkp))
       if os.access(self.__dumper_bkp2, os.F_OK): os.rmdir(self.__dumper_bkp2)
       return self.getCurrentDumper()


# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision