/[escript]/trunk-mpi-branch/escript/py_src/util.py
ViewVC logotype

Diff of /trunk-mpi-branch/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/py_src/util.py revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC trunk/escript/py_src/util.py revision 492 by gross, Fri Feb 3 02:07:24 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    #
3  ## @file util.py  #      COPYRIGHT ACcESS 2004 -  All Rights Reserved
4    #
5    #   This software is the property of ACcESS.  No part of this code
6    #   may be copied in any form or by any means without the expressed written
7    #   consent of ACcESS.  Copying, use or modification of this software
8    #   by any unauthorised person is illegal unless that
9    #   person has a software license agreement with ACcESS.
10    #
11    
12  """  """
13  @brief Utility functions for escript  Utility functions for escript
14    
15    @remark:  This module is under construction and is still tested!!!
16    
17    @var __author__: name of author
18    @var __licence__: licence agreement
19    @var __url__: url entry point on documentation
20    @var __version__: version
21    @var __date__: date of the version
22  """  """
23                                                                                                                                                                                                        
24    __author__="Lutz Gross, l.gross@uq.edu.au"
25    __licence__="contact: esys@access.uq.edu.au"
26    __url__="http://www.iservo.edu.au/esys/escript"
27    __version__="$Revision$"
28    __date__="$Date$"
29    
30    
31    import math
32  import numarray  import numarray
33    import numarray.linear_algebra
34  import escript  import escript
35    import os
36    
37    # missing tests:
38    
39    # def pokeShape(arg):
40    # def pokeDim(arg):
41    # def commonShape(arg0,arg1):
42    # def commonDim(*args):
43    # def testForZero(arg):
44    # def matchType(arg0=0.,arg1=0.):
45    # def matchShape(arg0,arg1):
46    
47    # def reorderComponents(arg,index):
48    
49  #  #
50  #   escript constants (have to be consistent witj utilC.h  # slicing: get
51  #  #          set
 UNKNOWN=-1  
 EPSILON=1.e-15  
 Pi=numarray.pi  
 # some solver options:  
 NO_REORDERING=0  
 MINIMUM_FILL_IN=1  
 NESTED_DISSECTION=2  
 # solver methods  
 DEFAULT_METHOD=0  
 DIRECT=1  
 CHOLEVSKY=2  
 PCG=3  
 CR=4  
 CGS=5  
 BICGSTAB=6  
 SSOR=7  
 ILU0=8  
 ILUT=9  
 JACOBI=10  
 GMRES=11  
 PRES20=12  
   
 METHOD_KEY="method"  
 SYMMETRY_KEY="symmetric"  
 TOLERANCE_KEY="tolerance"  
   
 # supported file formats:  
 VRML=1  
 PNG=2  
 JPEG=3  
 JPG=3  
 PS=4  
 OOGL=5  
 BMP=6  
 TIFF=7  
 OPENINVENTOR=8  
 RENDERMAN=9  
 PNM=10  
 #  
 # wrapper for various functions: if the argument has attribute the function name  
 # as an argument it calls the correspong methods. Otherwise the coresponsing numarray  
 # function is called.  
 #  
 # functions involving the underlying Domain:  
52  #  #
53  def grad(arg,where=None):  # and derivatives
54    
55    #=========================================================
56    #   some helpers:
57    #=========================================================
58    def saveVTK(filename,domain=None,**data):
59      """      """
60      @brief returns the spatial gradient of the Data object arg      writes a L{Data} objects into a files using the the VTK XML file format.
61    
62        Example:
63    
64      @param arg: Data object representing the function which gradient to be calculated.         tmp=Scalar(..)
65      @param where: FunctionSpace in which the gradient will be. If None Function(dom) where dom is the         v=Vector(..)
66                    domain of the Data object arg.         saveVTK("solution.xml",temperature=tmp,velovity=v)
67    
68        tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velovity"
69    
70        @param filename: file name of the output file
71        @type filename: C{str}
72        @param domain: domain of the L{Data} object. If not specified, the domain of the given L{Data} objects is used.
73        @type domain: L{escript.Domain}
74        @keyword <name>: writes the assigned value to the VTK file using <name> as identifier.
75        @type <name>: L{Data} object.
76        @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
77      """      """
78      if where==None:      if domain==None:
79         return arg.grad()         for i in data.keys():
80              if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()
81        if domain==None:
82            raise ValueError,"no domain detected."
83      else:      else:
84         return arg.grad(where)          domain.saveVTK(filename,data)
85    
86  def integrate(arg):  def saveDX(filename,domain=None,**data):
87      """      """
88      @brief return the integral if the function represented by Data object arg over its domain.      writes a L{Data} objects into a files using the the DX file format.
89    
90      @param arg      Example:
     """  
     return arg.integrate()  
91    
92  def interpolate(arg,where):         tmp=Scalar(..)
93           v=Vector(..)
94           saveDX("solution.dx",temperature=tmp,velovity=v)
95    
96        tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".
97    
98        @param filename: file name of the output file
99        @type filename: C{str}
100        @param domain: domain of the L{Data} object. If not specified, the domain of the given L{Data} objects is used.
101        @type domain: L{escript.Domain}
102        @keyword <name>: writes the assigned value to the DX file using <name> as identifier. The identifier can be used to select the data set when data are imported into DX.
103        @type <name>: L{Data} object.
104        @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
105      """      """
106      @brief interpolates the function represented by Data object arg into the FunctionSpace where.      if domain==None:
107           for i in data.keys():
108              if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()
109        if domain==None:
110            raise ValueError,"no domain detected."
111        else:
112            domain.saveDX(filename,data)
113    
114      @param arg  def kronecker(d=3):
115      @param where     """
116       return the kronecker S{delta}-symbol
117    
118       @param d: dimension or an object that has the C{getDim} method defining the dimension
119       @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
120       @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
121       @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2.
122       """
123       return identityTensor(d)
124    
125    def identity(shape=()):
126       """
127       return the shape x shape identity tensor
128    
129       @param shape: input shape for the identity tensor
130       @type shape: C{tuple} of C{int}
131       @return: array of shape shape x shape with M{u[i,k]=1} for M{i=k} and M{u[i,k]=0} otherwise for len(shape)=1 and, for len(shape)=2: M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise.
132       @rtype: L{numarray.NumArray} of rank 1, rankk 2 or rank 4.
133       @raise ValueError: if len(shape)>2.
134       """
135       if len(shape)>0:
136          out=numarray.zeros(shape+shape,numarray.Float)
137          if len(shape)==1:
138              for i0 in range(shape[0]):
139                 out[i0,i0]=1.
140          elif len(shape)==2:
141              for i0 in range(shape[0]):
142                 for i1 in range(shape[1]):
143                    out[i0,i1,i0,i1]=1.
144          else:
145              raise ValueError,"identity: length of shape is restricted to 2."
146       else:
147          out=1.
148       return out
149    
150    def identityTensor(d=3):
151       """
152       return the dxd identity matrix
153    
154       @param d: dimension or an object that has the C{getDim} method defining the dimension
155       @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
156       @return: the object u of rank 2 with M{u[i,j]=1} for M{i=j} and M{u[i,j]=0} otherwise
157       @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 2
158       """
159       if isinstance(d,escript.FunctionSpace):
160           return escript.Data(identity((d.getDim(),)),d)
161       elif isinstance(d,escript.Domain):
162           return identity((d.getDim(),))
163       else:
164           return identity((d,))
165    
166    def identityTensor4(d=3):
167       """
168       return the dxdxdxd identity tensor
169    
170       @param d: dimension or an object that has the C{getDim} method defining the dimension
171       @type d: C{int} or any object with a C{getDim} method
172       @return: the object u of rank 4 with M{u[i,j,k,l]=1} for M{i=k and j=l} and M{u[i,j,k,l]=0} otherwise
173       @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 4.
174       """
175       if isinstance(d,escript.FunctionSpace):
176           return escript.Data(identity((d.getDim(),d.getDim())),d)
177       elif isinstance(d,escript.Domain):
178           return identity((d.getDim(),d.getDim()))
179       else:
180           return identity((d,d))
181    
182    def unitVector(i=0,d=3):
183       """
184       return a unit vector u of dimension d with nonzero index i:
185    
186       @param i: index
187       @type i: C{int}
188       @param d: dimension or an object that has the C{getDim} method defining the dimension
189       @type d: C{int}, L{escript.Domain} or L{escript.FunctionSpace}
190       @return: the object u of rank 1 with M{u[j]=1} for M{j=i} and M{u[i]=0} otherwise
191       @rtype d: L{numarray.NumArray} or L{escript.Data} of rank 1
192       """
193       return kronecker(d)[i]
194    
195    #=========================================================================
196    #   global reduction operations (these functions have no symbolic version)
197    #=========================================================================
198    def Lsup(arg):
199      """      """
200      return arg.interpolate(where)      returns the Lsup-norm of argument arg. This is the maximum absolute value over all data points.
201        This function is equivalent to sup(abs(arg)).
202    
203  # functions returning Data objects:      @param arg: argument
204        @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
205        @return: maximum value of the absolute value of arg over all components and all data points
206        @rtype: C{float}
207        @raise TypeError: if type of arg cannot be processed
208        """
209        if isinstance(arg,numarray.NumArray):
210            return sup(abs(arg))
211        elif isinstance(arg,escript.Data):
212            return arg._Lsup()
213        elif isinstance(arg,float):
214            return abs(arg)
215        elif isinstance(arg,int):
216            return abs(float(arg))
217        else:
218          raise TypeError,"Lsup: Unknown argument type."
219    
220  def transpose(arg,axis=None):  def sup(arg):
221      """      """
222      @brief returns the transpose of the Data object arg.      returns the maximum value over all data points.
223    
224      @param arg      @param arg: argument
225        @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
226        @return: maximum value of arg over all components and all data points
227        @rtype: C{float}
228        @raise TypeError: if type of arg cannot be processed
229      """      """
230      if isinstance(arg,escript.Data):      if isinstance(arg,numarray.NumArray):
231         if axis==None: axis=arg.getRank()/2          return arg.max()
232         return arg.transpose(axis)      elif isinstance(arg,escript.Data):
233            return arg._sup()
234        elif isinstance(arg,float):
235            return arg
236        elif isinstance(arg,int):
237            return float(arg)
238      else:      else:
239         if axis==None: axis=arg.rank/2        raise TypeError,"sup: Unknown argument type."
        return numarray.transpose(arg,axis=axis)  
240    
241  def trace(arg):  def inf(arg):
242      """      """
243      @brief      returns the maximum value over all data points.
244    
245      @param arg      @param arg: argument
246        @type arg: C{float}, C{int}, L{escript.Data}, L{numarray.NumArray}.
247        @return : minimum value of arg over all components and all data points
248        @rtype: C{float}
249        @raise TypeError: if type of arg cannot be processed
250      """      """
251      if isinstance(arg,escript.Data):      if isinstance(arg,numarray.NumArray):
252         return arg.trace()          return arg.min()
253        elif isinstance(arg,escript.Data):
254            return arg._inf()
255        elif isinstance(arg,float):
256            return arg
257        elif isinstance(arg,int):
258            return float(arg)
259      else:      else:
260         return numarray.trace(arg)        raise TypeError,"inf: Unknown argument type."
261    
262  def exp(arg):  
263    #=========================================================================
264    #   some little helpers
265    #=========================================================================
266    def pokeShape(arg):
267      """      """
268      @brief      identifies the shape of its argument
269    
270      @param arg      @param arg: a given object
271        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
272        @return: the shape of the argument
273        @rtype: C{tuple} of C{int}
274        @raise TypeError: if type of arg cannot be processed
275      """      """
276      if isinstance(arg,escript.Data):  
277         return arg.exp()      if isinstance(arg,numarray.NumArray):
278            return arg.shape
279        elif isinstance(arg,escript.Data):
280            return arg.getShape()
281        elif isinstance(arg,float):
282            return ()
283        elif isinstance(arg,int):
284            return ()
285        elif isinstance(arg,Symbol):
286            return arg.getShape()
287      else:      else:
288         return numarray.exp(arg)        raise TypeError,"pokeShape: cannot identify shape"
289    
290  def sqrt(arg):  def pokeDim(arg):
291      """      """
292      @brief      identifies the spatial dimension of its argument
293    
294      @param arg      @param arg: a given object
295        @type arg: any
296        @return: the spatial dimension of the argument, if available, or C{None}
297        @rtype: C{int} or C{None}
298      """      """
299    
300      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
301         return arg.sqrt()          return arg.getFunctionSpace().getDim()
302        elif isinstance(arg,Symbol):
303            return arg.getDim()
304      else:      else:
305         return numarray.sqrt(arg)          return None
306    
307  def sin(arg):  def commonShape(arg0,arg1):
308      """      """
309      @brief      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
310    
311      @param arg      @param arg0: an object with a shape (see L{pokeShape})
312        @param arg1: an object with a shape (see L{pokeShape})
313        @return: the shape of arg0 or arg1 such that the left port equals the shape of arg0 and the right end equals the shape of arg1.
314        @rtype: C{tuple} of C{int}
315        @raise ValueError: if no shape can be found.
316      """      """
317      if isinstance(arg,escript.Data):      sh0=pokeShape(arg0)
318         return arg.sin()      sh1=pokeShape(arg1)
319        if len(sh0)<len(sh1):
320           if not sh0==sh1[:len(sh0)]:
321                 raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
322           return sh1
323        elif len(sh0)>len(sh1):
324           if not sh1==sh0[:len(sh1)]:
325                 raise ValueError,"argument 1 cannot be extended to the shape of argument 0"
326           return sh0
327      else:      else:
328         return numarray.sin(arg)         if not sh0==sh1:
329                 raise ValueError,"argument 1 and argument 0 have not the same shape."
330           return sh0
331    
332  def tan(arg):  def commonDim(*args):
333      """      """
334      @brief      identifies, if possible, the spatial dimension across a set of objects which may or my not have a spatial dimension.
335    
336      @param arg      @param *args: given objects
337        @return: the spatial dimension of the objects with identifiable dimension (see L{pokeDim}). If none the objects has
338                 a spatial dimension C{None} is returned.
339        @rtype: C{int} or C{None}
340        @raise ValueError: if the objects with identifiable dimension don't have the same spatial dimension.
341      """      """
342      if isinstance(arg,escript.Data):      out=None
343         return arg.tan()      for a in args:
344      else:         d=pokeDim(a)
345         return numarray.tan(arg)         if not out==None:
346              if not (d==None or out==d):
347                 raise ValueError,"dimension of arguments don't match"
348           else:
349              out=d
350        return out
351    
352  def cos(arg):  def testForZero(arg):
353      """      """
354      @brief      test the argument for being identical to Zero
355    
356      @param arg      @param arg: a given object
357        @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
358        @return : True if the argument is identical to zero.
359        @rtype : C{bool}
360      """      """
361      if isinstance(arg,escript.Data):      if isinstance(arg,numarray.NumArray):
362         return arg.cos()         return not Lsup(arg)>0.
363        elif isinstance(arg,escript.Data):
364           return False
365        elif isinstance(arg,float):
366           return not Lsup(arg)>0.
367        elif isinstance(arg,int):
368           return not Lsup(arg)>0.
369        elif isinstance(arg,Symbol):
370           return False
371      else:      else:
372         return numarray.cos(arg)         return False
373    
374  def maxval(arg):  def matchType(arg0=0.,arg1=0.):
375      """      """
376      @brief      converting arg0 and arg1 both to the same type L{numarray.NumArray} or L{escript.Data} or, if one of arg0 or arg1 is of type L{Symbol}, the other one to be of type L{numarray.NumArray} or L{escript.Data}.
377    
378      @param arg      @param arg0: first argument
379        @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
380        @param arg1: second argument
381        @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
382        @return: a tuple representing arg0 and arg1 with the same type or with one of them being a L{Symbol}
383        @rtype: C{tuple} of two L{numarray.NumArray}, two L{escript.Data}, a C{Symbol} and one of the types L{numarray.NumArray} or L{escript.Data}.  
384        @raise TypeError: if type of arg0 or arg1 cannot be processed
385      """      """
386      if isinstance(arg,escript.Data):      if isinstance(arg0,numarray.NumArray):
387         return arg.maxval()         if isinstance(arg1,numarray.NumArray):
388              pass
389           elif isinstance(arg1,escript.Data):
390              arg0=escript.Data(arg0,arg1.getFunctionSpace())
391           elif isinstance(arg1,float):
392              arg1=numarray.array(arg1)
393           elif isinstance(arg1,int):
394              arg1=numarray.array(float(arg1))
395           elif isinstance(arg1,Symbol):
396              pass
397           else:
398              raise TypeError,"function: Unknown type of second argument."    
399        elif isinstance(arg0,escript.Data):
400           if isinstance(arg1,numarray.NumArray):
401              arg1=escript.Data(arg1,arg0.getFunctionSpace())
402           elif isinstance(arg1,escript.Data):
403              pass
404           elif isinstance(arg1,float):
405              arg1=escript.Data(arg1,(),arg0.getFunctionSpace())
406           elif isinstance(arg1,int):
407              arg1=escript.Data(float(arg1),(),arg0.getFunctionSpace())
408           elif isinstance(arg1,Symbol):
409              pass
410           else:
411              raise TypeError,"function: Unknown type of second argument."    
412        elif isinstance(arg0,Symbol):
413           if isinstance(arg1,numarray.NumArray):
414              pass
415           elif isinstance(arg1,escript.Data):
416              pass
417           elif isinstance(arg1,float):
418              arg1=numarray.array(arg1)
419           elif isinstance(arg1,int):
420              arg1=numarray.array(float(arg1))
421           elif isinstance(arg1,Symbol):
422              pass
423           else:
424              raise TypeError,"function: Unknown type of second argument."    
425        elif isinstance(arg0,float):
426           if isinstance(arg1,numarray.NumArray):
427              arg0=numarray.array(arg0)
428           elif isinstance(arg1,escript.Data):
429              arg0=escript.Data(arg0,arg1.getFunctionSpace())
430           elif isinstance(arg1,float):
431              arg0=numarray.array(arg0)
432              arg1=numarray.array(arg1)
433           elif isinstance(arg1,int):
434              arg0=numarray.array(arg0)
435              arg1=numarray.array(float(arg1))
436           elif isinstance(arg1,Symbol):
437              arg0=numarray.array(arg0)
438           else:
439              raise TypeError,"function: Unknown type of second argument."    
440        elif isinstance(arg0,int):
441           if isinstance(arg1,numarray.NumArray):
442              arg0=numarray.array(float(arg0))
443           elif isinstance(arg1,escript.Data):
444              arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
445           elif isinstance(arg1,float):
446              arg0=numarray.array(float(arg0))
447              arg1=numarray.array(arg1)
448           elif isinstance(arg1,int):
449              arg0=numarray.array(float(arg0))
450              arg1=numarray.array(float(arg1))
451           elif isinstance(arg1,Symbol):
452              arg0=numarray.array(float(arg0))
453           else:
454              raise TypeError,"function: Unknown type of second argument."    
455      else:      else:
456         return arg.max()        raise TypeError,"function: Unknown type of first argument."    
457    
458  def minval(arg):      return arg0,arg1
459    
460    def matchShape(arg0,arg1):
461      """      """
462      @brief      
463    
464        If shape is not given the shape "largest" shape of args is used.
465    
466      @param arg      @param args: a given ob
467        @type arg: typically L{numarray.NumArray},L{escript.Data},C{float}, C{int}
468        @return: True if the argument is identical to zero.
469        @rtype: C{list} of C{int}
470      """      """
471      if isinstance(arg,escript.Data):      sh=commonShape(arg0,arg1)
472         return arg.minval()      sh0=pokeShape(arg0)
473      else:      sh1=pokeShape(arg1)
474         return arg.max()      if len(sh0)<len(sh):
475           return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1
476        elif len(sh1)<len(sh):
477           return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))
478        else:
479           return arg0,arg1
480    #=========================================================
481    #   symbolic tool box starts here:
482    #=========================================================
483    class Symbol(object):
484       """
485       Symbol class.
486    
487       Symbol class objects provide the same functionality as L{numarray.NumArray} and L{escript.Data} objects
488       but they do not have a value and therefore cannot be plotted or visualize. The main purpose is the possibilty
489       calculate derivatives with respect to other Symbols used to define a Symbol.
490    
491       """
492       def __init__(self,shape=(),args=[],dim=None):
493           """
494           Creates an instance of a symbol of a given shape. The symbol may depending on a list of arguments args which may be
495           symbols or any other object.
496    
497           @param arg: the arguments of the symbol.
498           @type arg: C{list}
499           @param shape: the shape
500           @type shape: C{tuple} of C{int}
501           @param dim: spatial dimension of the symbol. If dim=C{None} the spatial dimension is undefined.  
502           @type dim: C{None} or C{int}  
503    
504           """
505           if len(shape)>4:
506               raise ValueError,"Symbol supports only tensors up to order 4"
507           self.__args=args
508           self.__shape=shape
509           self.__dim=dim
510    
511       def getArgument(self,i=None):
512           """
513           returns the i-th argument of the symbol
514    
515           @param i: index of the argument requested.
516           @type i: C{int} or C{None}
517           @raise IndexError: if the requested index does not exist
518           @return: the vlaue of the i-th argument or i is not specified the list of all arguments.
519           @rtype: a single object or a list of objects
520           """
521           if i==None:
522              return self.__args
523           else:
524              if i<0 or i>=len(self.__args):
525                 raise IndexError,"there are only %s arguments"%len(self.__args)
526              return self.__args[i]
527    
528       def getRank(self):
529           """
530           the rank of the symbol
531    
532           @return: the rank of the symbol. This is length of the shape
533           @rtype: C{int}
534           """
535           return len(self.getShape())
536    
537       def getShape(self):
538           """
539           the shape of the symbol.
540    
541           @return : the shape of the symbol.
542           @rtype: C{tuple} of C{int}
543           """
544           return self.__shape
545    
546       def getDim(self):
547           """
548           the spatial dimension
549    
550           @return : the spatial dimension
551           @rtype: C{int} if the dimension is defined. Otherwise C{None} is returned.
552           """
553           return self.__dim
554      
555       def __str__(self):
556           """
557           a string representation of the symbol.
558           @return: a string representation of the object
559           @rtype: C{str}
560           """
561           args=[]
562           for arg in self.getArgument():
563              args.append(str(arg))
564           try:
565               out=self.getMyCode(args,format="str")
566           except NotImplementedError:
567               out="<Symbol %s>"%id(self)
568           return out
569          
570       def getSubstitutedArguments(self,argvals):
571           """
572           substitutes symbols in the arguments of this object and returns the result as a list.
573    
574           @param argvals: L{Symbols} and their substitutes. The L{Symbol} u in the expression defining this object is replaced by argvals[u].
575           @type argvals: C{dict} with keywords of type L{Symbol}.
576           @rtype: C{list} of objects
577           @return: list of the object assigned to the arguments through substitution or for the arguments which are not L{Symbols} the value assigned to the argument at instantiation.
578           """
579           out=[]
580           for a in self.getArgument():
581              if isinstance(a,Symbol):
582                 out.append(a.substitute(argvals))
583              else:
584                 out.append(a)
585           return out
586    
587       def getDifferentiatedArguments(self,arg):
588           """
589           applifies differentials to the arguments of this object and returns the result as a list.
590    
591           @param arg: the derivative is calculated with respect to arg
592           @type arg: typically L{escript.Symbol} but can also be C{float}, L{escript.Data}, L{numarray.NumArray} depending the involved functions and data.
593           @rtype: C{list} of objects
594           @return: list of object obtained by calculating the derivatives of the argumenst with respct to arg
595           """
596           out=[]
597           for a in self.getArgument():
598              if isinstance(a,Symbol):
599                 out.append(a.substitute(argvals))
600              else:
601                  s=pokeShape(s)+arg.getShape()
602                  if len(s)>0:
603                     out.append(numarray.zeros(s),numarray.Float)
604                  else:
605                     out.append(a)
606           return out
607    
608       def isAppropriateValue(self,arg):
609          """
610          checks if the given argument arg can be used as a substitution of this object. The method checks
611          the shape of arg and, if the spatial dimension is defined, the spatial dimension of arg.    
612    
613          @param arg: a given object
614          @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
615          @return: True if arg is a suitbale object to be used for substitution. Otherwise False is returned.
616          @rtype: C{bool}
617          """
618          if isinstance(arg,numarray.NumArray):
619              return arg.shape==self.getShape()
620          elif isinstance(arg,escript.Data):
621              if self.getDim()==None:
622                  return arg.getShape()==self.getShape()
623              elif self.getDim()==arg.getFunctionSpace().getDim():
624                  return arg.getShape()==self.getShape()
625              else:
626                  return False
627          elif isinstance(arg,Symbol):
628              if self.getDim()==None:
629                  return arg.getShape()==self.getShape()
630              elif self.getDim()==arg.getDim():
631                  return arg.getShape()==self.getShape()
632              else:
633                  return False
634          elif isinstance(arg,float):
635              return ()==self.getShape()
636          elif isinstance(arg,int):
637              return ()==self.getShape()
638          else:
639             return False
640    
641       def getMyCode(self,argstrs,format="escript"):
642           """
643           returns a program code that can be used to evaluate the symbol.
644    
645           @param argstrs: gives for each argument a string representing the argument for the evaluation.
646           @type argstrs: C{list} of C{str}.
647           @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
648           @type format: C{str}
649           @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
650           @rtype: C{str}
651           @raise NotImplementedError: if no implementation for the given format is available
652           @note: This method has to be overwritten by subclasses.
653           """
654           raise NotImplementedError,"no code for %s representation available"%format
655    
656       def substitute(self,argvals):
657          """  
658          assigns new values to symbols in the definition of the symbol.
659          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
660    
661          @param argvals: new values assigned to symbols
662          @type argvals: C{dict} with keywords of type L{Symbol}.
663          @return: result of the substitution process. Operations are executed as much as possible.
664          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
665          @note: this method has to be overwritten by a particular L{Symbol}
666          @raise NotImplementedError: if no implementation for the given format is available
667          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
668          """
669          if argvals.has_key(self):
670             arg=argvals[self]
671             if self.isAppropriateValue(arg):
672                return arg
673             else:
674                raise TypeError,"Symbol: new value is not appropriate."
675          else:
676             raise NotImplementedError,"no substitution in %s avialable"%str(self)
677    
678       def diff(self,arg):
679           """
680           returns the derivative of the symbol with respect to L{Symbol} arg
681    
682           @param arg: the derivative is calculated with respect to arg
683           @type arg: typically L{escript.Symbol} but can also be C{float}, L{escript.Data}, L{numarray.NumArray} depending the involved functions and data.
684           @return: derivative with respect to C{arg}
685           @rtype: typically L{escript.Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
686           @note: this method is overwritten by a particular L{Symbol}
687           """
688           if arg==self:
689              return identity(self.getShape())
690           else:
691              s=self.getShape()+arg.getShape()
692              if len(s)>0:
693                 return numarray.zeros(s,numarray.Float)
694              else:
695                 return 0.
696    
697       def __neg__(self):
698           """
699           returns -self.
700    
701           @return:  a S{Symbol} representing the negative of the object
702           @rtype: L{DependendSymbol}
703           """
704           return self*(-1.)
705    
706       def __pos__(self):
707           """
708           returns +self.
709    
710           @return:  a S{Symbol} representing the positive of the object
711           @rtype: L{DependendSymbol}
712           """
713           return self*(1.)
714    
715       def __abs__(self):
716           """
717           returns a S{Symbol} representing the absolute value of the object.
718           """
719           return Abs_Symbol(self)
720    
721       def __add__(self,other):
722           """
723           add another object to this object
724    
725           @param other: object to be added to this object
726           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
727           @return:  a S{Symbol} representing the sum of this object and C{other}
728           @rtype: L{DependendSymbol}
729           """
730           return add(self,other)
731    
732       def __radd__(self,other):
733           """
734           add this object to another object
735    
736           @param other: object this object is added to
737           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
738           @return: a S{Symbol} representing the sum of C{other} and this object object
739           @rtype: L{DependendSymbol}
740           """
741           return add(other,self)
742    
743       def __sub__(self,other):
744           """
745           subtracts another object from this object
746    
747           @param other: object to be subtracted from this object
748           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
749           @return: a S{Symbol} representing the difference of C{other} and this object
750           @rtype: L{DependendSymbol}
751           """
752           return add(self,-other)
753    
754       def __rsub__(self,other):
755           """
756           subtracts this object from another object
757    
758           @param other: object this object is been subtracted from
759           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
760           @return: a S{Symbol} representing the difference of this object and C{other}.
761           @rtype: L{DependendSymbol}
762           """
763           return add(-self,other)
764    
765       def __mul__(self,other):
766           """
767           multiplies this object with other object
768    
769           @param other: object to be mutiplied by this object
770           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
771           @return: a S{Symbol} representing the product of the object and C{other}.
772           @rtype: L{DependendSymbol} or 0 if other is identical to zero.
773           """
774           return mult(self,other)
775    
776       def __rmul__(self,other):
777           """
778           multiplies this object with other object
779    
780           @param other: object this object is multiplied with
781           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
782           @return: a S{Symbol} representing the product of C{other} and the object.
783           @rtype: L{DependendSymbol} or 0 if other is identical to zero.
784           """
785           return mult(other,self)
786    
787       def __div__(self,other):
788           """
789           divides this object by other object
790    
791           @param other: object dividing this object
792           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
793           @return: a S{Symbol} representing the quotient of this object and C{other}
794           @rtype: L{DependendSymbol}
795           """
796           return quotient(self,other)
797    
798       def __rdiv__(self,other):
799           """
800           divides this object by other object
801    
802           @param other: object dividing this object
803           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
804           @return: a S{Symbol} representing the quotient of C{other} and this object
805           @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
806           """
807           return quotient(other,self)
808    
809       def __pow__(self,other):
810           """
811           raises this object to the power of other
812    
813           @param other: exponent
814           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
815           @return: a S{Symbol} representing the power of this object to C{other}
816           @rtype: L{DependendSymbol} or 1 if C{other} is identical to zero.
817           """
818           return power(self,other)
819    
820       def __rpow__(self,other):
821           """
822           raises an object to the power of this object
823    
824           @param other: basis
825           @type other: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
826           @return: a S{Symbol} representing the power of C{other} to this object
827           @rtype: L{DependendSymbol} or 0 if C{other} is identical to zero.
828           """
829           return power(other,self)
830    
831    class DependendSymbol(Symbol):
832       """
833       DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.
834       Two DependendSymbol are equal if they have the same shape, the same arguments and one of them has an unspecified spatial dimension or the spatial dimension is identical  
835      
836       Example:
837      
838       u1=Symbol(shape=(3,4),dim=2,args=[4.])
839       u2=Symbol(shape=(3,4),dim=2,args=[4.])
840       print u1==u2
841       False
842      
843          but  
844    
845       u1=DependendSymbol(shape=(3,4),dim=2,args=[4.])
846       u2=DependendSymbol(shape=(3,4),dim=2,args=[4.])
847       u3=DependendSymbol(shape=(2,),dim=2,args=[4.])  
848       print u1==u2, u1==u3
849       True False
850    
851       @note: DependendSymbol should be used as return value of functions with L{Symbol} arguments. This will allow the optimizer to remove redundant function calls.
852       """
853       def __eq__(self,other):
854          """
855          checks if other equals self
856    
857          @param other: any object
858          @return: True if other has the same class like self, and the shape, the spatial diemsnion and the arguments are equal.
859          @rtype: C{bool}
860          """
861          if isinstance(other,DependendSymbol):
862             if self.__class__==other.__class__:  
863                if self.getShape()==other.getShape():
864                  if self.getArgument()==other.getArgument():
865                     if self.getDim()==None or other.getDim()==None or self.getDim()==other.getDim():
866                        return True
867          return False
868    
869       def __ne__(self,other):
870          """
871          checks if other equals self
872    
873          @param other: any object
874          @return: Flase if other has the same class like self, and the shape, the spatial diemsnion and the arguments are equal.
875          @rtype: C{bool}
876          """
877          return not self==other
878    #=========================================================
879    #  Unary operations prserving the shape
880    #========================================================
881    def log10(arg):
882       """
883       returns base-10 logarithm of argument arg
884    
885       @param arg: argument
886       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
887       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
888       @raises TypeError: if the type of the argument is not expected.
889       """
890       if isinstance(arg,numarray.NumArray):
891          return numarray.log10(arg)
892       elif isinstance(arg,escript.Data):
893          return arg._log10()
894       elif isinstance(arg,float):
895          return math.log10(arg)
896       elif isinstance(arg,int):
897          return math.log10(float(arg))
898       elif isinstance(arg,Symbol):
899          return log(arg)/log(10.)
900       else:
901          raise TypeError,"log10: Unknown argument type."
902    
903    def wherePositive(arg):
904       """
905       returns mask of positive values of argument arg
906    
907       @param arg: argument
908       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
909       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
910       @raises TypeError: if the type of the argument is not expected.
911       """
912       if isinstance(arg,numarray.NumArray):
913          out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))*1.
914          if isinstance(out,float): out=numarray.array(out)
915          return out
916       elif isinstance(arg,escript.Data):
917          return arg._wherePositive()
918       elif isinstance(arg,float):
919          if arg>0:
920            return 1.
921          else:
922            return 0.
923       elif isinstance(arg,int):
924          if arg>0:
925            return 1.
926          else:
927            return 0.
928       elif isinstance(arg,Symbol):
929          return WherePositive_Symbol(arg)
930       else:
931          raise TypeError,"wherePositive: Unknown argument type."
932    
933    class WherePositive_Symbol(DependendSymbol):
934       """
935       L{Symbol} representing the result of the mask of positive values function
936       """
937       def __init__(self,arg):
938          """
939          initialization of wherePositive L{Symbol} with argument arg
940          @param arg: argument of function
941          @type arg: typically L{Symbol}.
942          """
943          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
944    
945       def getMyCode(self,argstrs,format="escript"):
946          """
947          returns a program code that can be used to evaluate the symbol.
948    
949          @param argstrs: gives for each argument a string representing the argument for the evaluation.
950          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
951          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
952          @type format: C{str}
953          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
954          @rtype: C{str}
955          @raise: NotImplementedError: if the requested format is not available
956          """
957          if isinstance(argstrs,list):
958              argstrs=argstrs[0]
959          if format=="escript" or format=="str"  or format=="text":
960             return "wherePositive(%s)"%argstrs
961          else:
962             raise NotImplementedError,"WherePositive_Symbol does not provide program code for format %s."%format
963    
964       def substitute(self,argvals):
965          """
966          assigns new values to symbols in the definition of the symbol.
967          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
968    
969          @param argvals: new values assigned to symbols
970          @type argvals: C{dict} with keywords of type L{Symbol}.
971          @return: result of the substitution process. Operations are executed as much as possible.
972          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
973          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
974          """
975          if argvals.has_key(self):
976             arg=argvals[self]
977             if self.isAppropriateValue(arg):
978                return arg
979             else:
980                raise TypeError,"%s: new value is not appropriate."%str(self)
981          else:
982             arg=self.getSubstitutedArguments(argvals)[0]
983             return wherePositive(arg)
984    
985    def whereNegative(arg):
986       """
987       returns mask of positive values of argument arg
988    
989       @param arg: argument
990       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
991       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
992       @raises TypeError: if the type of the argument is not expected.
993       """
994       if isinstance(arg,numarray.NumArray):
995          out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))*1.
996          if isinstance(out,float): out=numarray.array(out)
997          return out
998       elif isinstance(arg,escript.Data):
999          return arg._whereNegative()
1000       elif isinstance(arg,float):
1001          if arg<0:
1002            return 1.
1003          else:
1004            return 0.
1005       elif isinstance(arg,int):
1006          if arg<0:
1007            return 1.
1008          else:
1009            return 0.
1010       elif isinstance(arg,Symbol):
1011          return WhereNegative_Symbol(arg)
1012       else:
1013          raise TypeError,"whereNegative: Unknown argument type."
1014    
1015    class WhereNegative_Symbol(DependendSymbol):
1016       """
1017       L{Symbol} representing the result of the mask of positive values function
1018       """
1019       def __init__(self,arg):
1020          """
1021          initialization of whereNegative L{Symbol} with argument arg
1022          @param arg: argument of function
1023          @type arg: typically L{Symbol}.
1024          """
1025          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1026    
1027       def getMyCode(self,argstrs,format="escript"):
1028          """
1029          returns a program code that can be used to evaluate the symbol.
1030    
1031          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1032          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1033          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1034          @type format: C{str}
1035          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1036          @rtype: C{str}
1037          @raise: NotImplementedError: if the requested format is not available
1038          """
1039          if isinstance(argstrs,list):
1040              argstrs=argstrs[0]
1041          if format=="escript" or format=="str"  or format=="text":
1042             return "whereNegative(%s)"%argstrs
1043          else:
1044             raise NotImplementedError,"WhereNegative_Symbol does not provide program code for format %s."%format
1045    
1046       def substitute(self,argvals):
1047          """
1048          assigns new values to symbols in the definition of the symbol.
1049          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1050    
1051          @param argvals: new values assigned to symbols
1052          @type argvals: C{dict} with keywords of type L{Symbol}.
1053          @return: result of the substitution process. Operations are executed as much as possible.
1054          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1055          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1056          """
1057          if argvals.has_key(self):
1058             arg=argvals[self]
1059             if self.isAppropriateValue(arg):
1060                return arg
1061             else:
1062                raise TypeError,"%s: new value is not appropriate."%str(self)
1063          else:
1064             arg=self.getSubstitutedArguments(argvals)[0]
1065             return whereNegative(arg)
1066    
1067    def whereNonNegative(arg):
1068       """
1069       returns mask of non-negative values of argument arg
1070    
1071       @param arg: argument
1072       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1073       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1074       @raises TypeError: if the type of the argument is not expected.
1075       """
1076       if isinstance(arg,numarray.NumArray):
1077          out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.
1078          if isinstance(out,float): out=numarray.array(out)
1079          return out
1080       elif isinstance(arg,escript.Data):
1081          return arg._whereNonNegative()
1082       elif isinstance(arg,float):
1083          if arg<0:
1084            return 0.
1085          else:
1086            return 1.
1087       elif isinstance(arg,int):
1088          if arg<0:
1089            return 0.
1090          else:
1091            return 1.
1092       elif isinstance(arg,Symbol):
1093          return 1.-whereNegative(arg)
1094       else:
1095          raise TypeError,"whereNonNegative: Unknown argument type."
1096    
1097    def whereNonPositive(arg):
1098       """
1099       returns mask of non-positive values of argument arg
1100    
1101       @param arg: argument
1102       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1103       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1104       @raises TypeError: if the type of the argument is not expected.
1105       """
1106       if isinstance(arg,numarray.NumArray):
1107          out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.
1108          if isinstance(out,float): out=numarray.array(out)
1109          return out
1110       elif isinstance(arg,escript.Data):
1111          return arg._whereNonPositive()
1112       elif isinstance(arg,float):
1113          if arg>0:
1114            return 0.
1115          else:
1116            return 1.
1117       elif isinstance(arg,int):
1118          if arg>0:
1119            return 0.
1120          else:
1121            return 1.
1122       elif isinstance(arg,Symbol):
1123          return 1.-wherePositive(arg)
1124       else:
1125          raise TypeError,"whereNonPositive: Unknown argument type."
1126    
1127    def whereZero(arg,tol=0.):
1128       """
1129       returns mask of zero entries of argument arg
1130    
1131       @param arg: argument
1132       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1133       @param tol: tolerance. values with absolute value less then tol are accepted as zero.
1134       @type tol: C{float}
1135       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1136       @raises TypeError: if the type of the argument is not expected.
1137       """
1138       if isinstance(arg,numarray.NumArray):
1139          out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.
1140          if isinstance(out,float): out=numarray.array(out)
1141          return out
1142       elif isinstance(arg,escript.Data):
1143          if tol>0.:
1144             return whereNegative(abs(arg)-tol)
1145          else:
1146             return arg._whereZero()
1147       elif isinstance(arg,float):
1148          if abs(arg)<=tol:
1149            return 1.
1150          else:
1151            return 0.
1152       elif isinstance(arg,int):
1153          if abs(float(arg))<=tol:
1154            return 1.
1155          else:
1156            return 0.
1157       elif isinstance(arg,Symbol):
1158          return WhereZero_Symbol(arg,tol)
1159       else:
1160          raise TypeError,"whereZero: Unknown argument type."
1161    
1162    class WhereZero_Symbol(DependendSymbol):
1163       """
1164       L{Symbol} representing the result of the mask of zero entries function
1165       """
1166       def __init__(self,arg,tol=0.):
1167          """
1168          initialization of whereZero L{Symbol} with argument arg
1169          @param arg: argument of function
1170          @type arg: typically L{Symbol}.
1171          """
1172          DependendSymbol.__init__(self,args=[arg,tol],shape=arg.getShape(),dim=arg.getDim())
1173    
1174       def getMyCode(self,argstrs,format="escript"):
1175          """
1176          returns a program code that can be used to evaluate the symbol.
1177    
1178          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1179          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1180          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1181          @type format: C{str}
1182          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1183          @rtype: C{str}
1184          @raise: NotImplementedError: if the requested format is not available
1185          """
1186          if format=="escript" or format=="str"  or format=="text":
1187             return "whereZero(%s,tol=%s)"%(argstrs[0],argstrs[1])
1188          else:
1189             raise NotImplementedError,"WhereZero_Symbol does not provide program code for format %s."%format
1190    
1191       def substitute(self,argvals):
1192          """
1193          assigns new values to symbols in the definition of the symbol.
1194          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1195    
1196          @param argvals: new values assigned to symbols
1197          @type argvals: C{dict} with keywords of type L{Symbol}.
1198          @return: result of the substitution process. Operations are executed as much as possible.
1199          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1200          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1201          """
1202          if argvals.has_key(self):
1203             arg=argvals[self]
1204             if self.isAppropriateValue(arg):
1205                return arg
1206             else:
1207                raise TypeError,"%s: new value is not appropriate."%str(self)
1208          else:
1209             arg=self.getSubstitutedArguments(argvals)
1210             return whereZero(arg[0],arg[1])
1211    
1212    def whereNonZero(arg,tol=0.):
1213       """
1214       returns mask of values different from zero of argument arg
1215    
1216       @param arg: argument
1217       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1218       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1219       @raises TypeError: if the type of the argument is not expected.
1220       """
1221       if isinstance(arg,numarray.NumArray):
1222          out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.
1223          if isinstance(out,float): out=numarray.array(out)
1224          return out
1225       elif isinstance(arg,escript.Data):
1226          if tol>0.:
1227             return 1.-whereZero(arg,tol)
1228          else:
1229             return arg._whereNonZero()
1230       elif isinstance(arg,float):
1231          if abs(arg)>tol:
1232            return 1.
1233          else:
1234            return 0.
1235       elif isinstance(arg,int):
1236          if abs(float(arg))>tol:
1237            return 1.
1238          else:
1239            return 0.
1240       elif isinstance(arg,Symbol):
1241          return 1.-whereZero(arg,tol)
1242       else:
1243          raise TypeError,"whereNonZero: Unknown argument type."
1244    
1245    def sin(arg):
1246       """
1247       returns sine of argument arg
1248    
1249       @param arg: argument
1250       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1251       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1252       @raises TypeError: if the type of the argument is not expected.
1253       """
1254       if isinstance(arg,numarray.NumArray):
1255          return numarray.sin(arg)
1256       elif isinstance(arg,escript.Data):
1257          return arg._sin()
1258       elif isinstance(arg,float):
1259          return math.sin(arg)
1260       elif isinstance(arg,int):
1261          return math.sin(arg)
1262       elif isinstance(arg,Symbol):
1263          return Sin_Symbol(arg)
1264       else:
1265          raise TypeError,"sin: Unknown argument type."
1266    
1267    class Sin_Symbol(DependendSymbol):
1268       """
1269       L{Symbol} representing the result of the sine function
1270       """
1271       def __init__(self,arg):
1272          """
1273          initialization of sin L{Symbol} with argument arg
1274          @param arg: argument of function
1275          @type arg: typically L{Symbol}.
1276          """
1277          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1278    
1279       def getMyCode(self,argstrs,format="escript"):
1280          """
1281          returns a program code that can be used to evaluate the symbol.
1282    
1283          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1284          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1285          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1286          @type format: C{str}
1287          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1288          @rtype: C{str}
1289          @raise: NotImplementedError: if the requested format is not available
1290          """
1291          if isinstance(argstrs,list):
1292              argstrs=argstrs[0]
1293          if format=="escript" or format=="str"  or format=="text":
1294             return "sin(%s)"%argstrs
1295          else:
1296             raise NotImplementedError,"Sin_Symbol does not provide program code for format %s."%format
1297    
1298       def substitute(self,argvals):
1299          """
1300          assigns new values to symbols in the definition of the symbol.
1301          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1302    
1303          @param argvals: new values assigned to symbols
1304          @type argvals: C{dict} with keywords of type L{Symbol}.
1305          @return: result of the substitution process. Operations are executed as much as possible.
1306          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1307          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1308          """
1309          if argvals.has_key(self):
1310             arg=argvals[self]
1311             if self.isAppropriateValue(arg):
1312                return arg
1313             else:
1314                raise TypeError,"%s: new value is not appropriate."%str(self)
1315          else:
1316             arg=self.getSubstitutedArguments(argvals)[0]
1317             return sin(arg)
1318    
1319       def diff(self,arg):
1320          """
1321          differential of this object
1322    
1323          @param arg: the derivative is calculated with respect to arg
1324          @type arg: L{escript.Symbol}
1325          @return: derivative with respect to C{arg}
1326          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1327          """
1328          if arg==self:
1329             return identity(self.getShape())
1330          else:
1331             myarg=self.getArgument()[0]
1332             val=matchShape(cos(myarg),self.getDifferentiatedArguments(arg)[0])
1333             return val[0]*val[1]
1334    
1335    def cos(arg):
1336       """
1337       returns cosine of argument arg
1338    
1339       @param arg: argument
1340       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1341       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1342       @raises TypeError: if the type of the argument is not expected.
1343       """
1344       if isinstance(arg,numarray.NumArray):
1345          return numarray.cos(arg)
1346       elif isinstance(arg,escript.Data):
1347          return arg._cos()
1348       elif isinstance(arg,float):
1349          return math.cos(arg)
1350       elif isinstance(arg,int):
1351          return math.cos(arg)
1352       elif isinstance(arg,Symbol):
1353          return Cos_Symbol(arg)
1354       else:
1355          raise TypeError,"cos: Unknown argument type."
1356    
1357    class Cos_Symbol(DependendSymbol):
1358       """
1359       L{Symbol} representing the result of the cosine function
1360       """
1361       def __init__(self,arg):
1362          """
1363          initialization of cos L{Symbol} with argument arg
1364          @param arg: argument of function
1365          @type arg: typically L{Symbol}.
1366          """
1367          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1368    
1369       def getMyCode(self,argstrs,format="escript"):
1370          """
1371          returns a program code that can be used to evaluate the symbol.
1372    
1373          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1374          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1375          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1376          @type format: C{str}
1377          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1378          @rtype: C{str}
1379          @raise: NotImplementedError: if the requested format is not available
1380          """
1381          if isinstance(argstrs,list):
1382              argstrs=argstrs[0]
1383          if format=="escript" or format=="str"  or format=="text":
1384             return "cos(%s)"%argstrs
1385          else:
1386             raise NotImplementedError,"Cos_Symbol does not provide program code for format %s."%format
1387    
1388       def substitute(self,argvals):
1389          """
1390          assigns new values to symbols in the definition of the symbol.
1391          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1392    
1393          @param argvals: new values assigned to symbols
1394          @type argvals: C{dict} with keywords of type L{Symbol}.
1395          @return: result of the substitution process. Operations are executed as much as possible.
1396          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1397          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1398          """
1399          if argvals.has_key(self):
1400             arg=argvals[self]
1401             if self.isAppropriateValue(arg):
1402                return arg
1403             else:
1404                raise TypeError,"%s: new value is not appropriate."%str(self)
1405          else:
1406             arg=self.getSubstitutedArguments(argvals)[0]
1407             return cos(arg)
1408    
1409       def diff(self,arg):
1410          """
1411          differential of this object
1412    
1413          @param arg: the derivative is calculated with respect to arg
1414          @type arg: L{escript.Symbol}
1415          @return: derivative with respect to C{arg}
1416          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1417          """
1418          if arg==self:
1419             return identity(self.getShape())
1420          else:
1421             myarg=self.getArgument()[0]
1422             val=matchShape(-sin(myarg),self.getDifferentiatedArguments(arg)[0])
1423             return val[0]*val[1]
1424    
1425    def tan(arg):
1426       """
1427       returns tangent of argument arg
1428    
1429       @param arg: argument
1430       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1431       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1432       @raises TypeError: if the type of the argument is not expected.
1433       """
1434       if isinstance(arg,numarray.NumArray):
1435          return numarray.tan(arg)
1436       elif isinstance(arg,escript.Data):
1437          return arg._tan()
1438       elif isinstance(arg,float):
1439          return math.tan(arg)
1440       elif isinstance(arg,int):
1441          return math.tan(arg)
1442       elif isinstance(arg,Symbol):
1443          return Tan_Symbol(arg)
1444       else:
1445          raise TypeError,"tan: Unknown argument type."
1446    
1447    class Tan_Symbol(DependendSymbol):
1448       """
1449       L{Symbol} representing the result of the tangent function
1450       """
1451       def __init__(self,arg):
1452          """
1453          initialization of tan L{Symbol} with argument arg
1454          @param arg: argument of function
1455          @type arg: typically L{Symbol}.
1456          """
1457          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1458    
1459       def getMyCode(self,argstrs,format="escript"):
1460          """
1461          returns a program code that can be used to evaluate the symbol.
1462    
1463          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1464          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1465          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1466          @type format: C{str}
1467          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1468          @rtype: C{str}
1469          @raise: NotImplementedError: if the requested format is not available
1470          """
1471          if isinstance(argstrs,list):
1472              argstrs=argstrs[0]
1473          if format=="escript" or format=="str"  or format=="text":
1474             return "tan(%s)"%argstrs
1475          else:
1476             raise NotImplementedError,"Tan_Symbol does not provide program code for format %s."%format
1477    
1478       def substitute(self,argvals):
1479          """
1480          assigns new values to symbols in the definition of the symbol.
1481          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1482    
1483          @param argvals: new values assigned to symbols
1484          @type argvals: C{dict} with keywords of type L{Symbol}.
1485          @return: result of the substitution process. Operations are executed as much as possible.
1486          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1487          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1488          """
1489          if argvals.has_key(self):
1490             arg=argvals[self]
1491             if self.isAppropriateValue(arg):
1492                return arg
1493             else:
1494                raise TypeError,"%s: new value is not appropriate."%str(self)
1495          else:
1496             arg=self.getSubstitutedArguments(argvals)[0]
1497             return tan(arg)
1498    
1499       def diff(self,arg):
1500          """
1501          differential of this object
1502    
1503          @param arg: the derivative is calculated with respect to arg
1504          @type arg: L{escript.Symbol}
1505          @return: derivative with respect to C{arg}
1506          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1507          """
1508          if arg==self:
1509             return identity(self.getShape())
1510          else:
1511             myarg=self.getArgument()[0]
1512             val=matchShape(1./cos(myarg)**2,self.getDifferentiatedArguments(arg)[0])
1513             return val[0]*val[1]
1514    
1515    def asin(arg):
1516       """
1517       returns inverse sine of argument arg
1518    
1519       @param arg: argument
1520       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1521       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1522       @raises TypeError: if the type of the argument is not expected.
1523       """
1524       if isinstance(arg,numarray.NumArray):
1525          return numarray.arcsin(arg)
1526       elif isinstance(arg,escript.Data):
1527          return arg._asin()
1528       elif isinstance(arg,float):
1529          return math.asin(arg)
1530       elif isinstance(arg,int):
1531          return math.asin(arg)
1532       elif isinstance(arg,Symbol):
1533          return Asin_Symbol(arg)
1534       else:
1535          raise TypeError,"asin: Unknown argument type."
1536    
1537    class Asin_Symbol(DependendSymbol):
1538       """
1539       L{Symbol} representing the result of the inverse sine function
1540       """
1541       def __init__(self,arg):
1542          """
1543          initialization of asin L{Symbol} with argument arg
1544          @param arg: argument of function
1545          @type arg: typically L{Symbol}.
1546          """
1547          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1548    
1549       def getMyCode(self,argstrs,format="escript"):
1550          """
1551          returns a program code that can be used to evaluate the symbol.
1552    
1553          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1554          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1555          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1556          @type format: C{str}
1557          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1558          @rtype: C{str}
1559          @raise: NotImplementedError: if the requested format is not available
1560          """
1561          if isinstance(argstrs,list):
1562              argstrs=argstrs[0]
1563          if format=="escript" or format=="str"  or format=="text":
1564             return "asin(%s)"%argstrs
1565          else:
1566             raise NotImplementedError,"Asin_Symbol does not provide program code for format %s."%format
1567    
1568       def substitute(self,argvals):
1569          """
1570          assigns new values to symbols in the definition of the symbol.
1571          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1572    
1573          @param argvals: new values assigned to symbols
1574          @type argvals: C{dict} with keywords of type L{Symbol}.
1575          @return: result of the substitution process. Operations are executed as much as possible.
1576          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1577          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1578          """
1579          if argvals.has_key(self):
1580             arg=argvals[self]
1581             if self.isAppropriateValue(arg):
1582                return arg
1583             else:
1584                raise TypeError,"%s: new value is not appropriate."%str(self)
1585          else:
1586             arg=self.getSubstitutedArguments(argvals)[0]
1587             return asin(arg)
1588    
1589       def diff(self,arg):
1590          """
1591          differential of this object
1592    
1593          @param arg: the derivative is calculated with respect to arg
1594          @type arg: L{escript.Symbol}
1595          @return: derivative with respect to C{arg}
1596          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1597          """
1598          if arg==self:
1599             return identity(self.getShape())
1600          else:
1601             myarg=self.getArgument()[0]
1602             val=matchShape(1./sqrt(1.-myarg**2),self.getDifferentiatedArguments(arg)[0])
1603             return val[0]*val[1]
1604    
1605    def acos(arg):
1606       """
1607       returns inverse cosine of argument arg
1608    
1609       @param arg: argument
1610       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1611       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1612       @raises TypeError: if the type of the argument is not expected.
1613       """
1614       if isinstance(arg,numarray.NumArray):
1615          return numarray.arccos(arg)
1616       elif isinstance(arg,escript.Data):
1617          return arg._acos()
1618       elif isinstance(arg,float):
1619          return math.acos(arg)
1620       elif isinstance(arg,int):
1621          return math.acos(arg)
1622       elif isinstance(arg,Symbol):
1623          return Acos_Symbol(arg)
1624       else:
1625          raise TypeError,"acos: Unknown argument type."
1626    
1627    class Acos_Symbol(DependendSymbol):
1628       """
1629       L{Symbol} representing the result of the inverse cosine function
1630       """
1631       def __init__(self,arg):
1632          """
1633          initialization of acos L{Symbol} with argument arg
1634          @param arg: argument of function
1635          @type arg: typically L{Symbol}.
1636          """
1637          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1638    
1639       def getMyCode(self,argstrs,format="escript"):
1640          """
1641          returns a program code that can be used to evaluate the symbol.
1642    
1643          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1644          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1645          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1646          @type format: C{str}
1647          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1648          @rtype: C{str}
1649          @raise: NotImplementedError: if the requested format is not available
1650          """
1651          if isinstance(argstrs,list):
1652              argstrs=argstrs[0]
1653          if format=="escript" or format=="str"  or format=="text":
1654             return "acos(%s)"%argstrs
1655          else:
1656             raise NotImplementedError,"Acos_Symbol does not provide program code for format %s."%format
1657    
1658       def substitute(self,argvals):
1659          """
1660          assigns new values to symbols in the definition of the symbol.
1661          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1662    
1663          @param argvals: new values assigned to symbols
1664          @type argvals: C{dict} with keywords of type L{Symbol}.
1665          @return: result of the substitution process. Operations are executed as much as possible.
1666          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1667          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1668          """
1669          if argvals.has_key(self):
1670             arg=argvals[self]
1671             if self.isAppropriateValue(arg):
1672                return arg
1673             else:
1674                raise TypeError,"%s: new value is not appropriate."%str(self)
1675          else:
1676             arg=self.getSubstitutedArguments(argvals)[0]
1677             return acos(arg)
1678    
1679       def diff(self,arg):
1680          """
1681          differential of this object
1682    
1683          @param arg: the derivative is calculated with respect to arg
1684          @type arg: L{escript.Symbol}
1685          @return: derivative with respect to C{arg}
1686          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1687          """
1688          if arg==self:
1689             return identity(self.getShape())
1690          else:
1691             myarg=self.getArgument()[0]
1692             val=matchShape(-1./sqrt(1.-myarg**2),self.getDifferentiatedArguments(arg)[0])
1693             return val[0]*val[1]
1694    
1695    def atan(arg):
1696       """
1697       returns inverse tangent of argument arg
1698    
1699       @param arg: argument
1700       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1701       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1702       @raises TypeError: if the type of the argument is not expected.
1703       """
1704       if isinstance(arg,numarray.NumArray):
1705          return numarray.arctan(arg)
1706       elif isinstance(arg,escript.Data):
1707          return arg._atan()
1708       elif isinstance(arg,float):
1709          return math.atan(arg)
1710       elif isinstance(arg,int):
1711          return math.atan(arg)
1712       elif isinstance(arg,Symbol):
1713          return Atan_Symbol(arg)
1714       else:
1715          raise TypeError,"atan: Unknown argument type."
1716    
1717    class Atan_Symbol(DependendSymbol):
1718       """
1719       L{Symbol} representing the result of the inverse tangent function
1720       """
1721       def __init__(self,arg):
1722          """
1723          initialization of atan L{Symbol} with argument arg
1724          @param arg: argument of function
1725          @type arg: typically L{Symbol}.
1726          """
1727          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1728    
1729       def getMyCode(self,argstrs,format="escript"):
1730          """
1731          returns a program code that can be used to evaluate the symbol.
1732    
1733          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1734          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1735          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1736          @type format: C{str}
1737          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1738          @rtype: C{str}
1739          @raise: NotImplementedError: if the requested format is not available
1740          """
1741          if isinstance(argstrs,list):
1742              argstrs=argstrs[0]
1743          if format=="escript" or format=="str"  or format=="text":
1744             return "atan(%s)"%argstrs
1745          else:
1746             raise NotImplementedError,"Atan_Symbol does not provide program code for format %s."%format
1747    
1748       def substitute(self,argvals):
1749          """
1750          assigns new values to symbols in the definition of the symbol.
1751          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1752    
1753          @param argvals: new values assigned to symbols
1754          @type argvals: C{dict} with keywords of type L{Symbol}.
1755          @return: result of the substitution process. Operations are executed as much as possible.
1756          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1757          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1758          """
1759          if argvals.has_key(self):
1760             arg=argvals[self]
1761             if self.isAppropriateValue(arg):
1762                return arg
1763             else:
1764                raise TypeError,"%s: new value is not appropriate."%str(self)
1765          else:
1766             arg=self.getSubstitutedArguments(argvals)[0]
1767             return atan(arg)
1768    
1769       def diff(self,arg):
1770          """
1771          differential of this object
1772    
1773          @param arg: the derivative is calculated with respect to arg
1774          @type arg: L{escript.Symbol}
1775          @return: derivative with respect to C{arg}
1776          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1777          """
1778          if arg==self:
1779             return identity(self.getShape())
1780          else:
1781             myarg=self.getArgument()[0]
1782             val=matchShape(1./(1+myarg**2),self.getDifferentiatedArguments(arg)[0])
1783             return val[0]*val[1]
1784    
1785    def sinh(arg):
1786       """
1787       returns hyperbolic sine of argument arg
1788    
1789       @param arg: argument
1790       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1791       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1792       @raises TypeError: if the type of the argument is not expected.
1793       """
1794       if isinstance(arg,numarray.NumArray):
1795          return numarray.sinh(arg)
1796       elif isinstance(arg,escript.Data):
1797          return arg._sinh()
1798       elif isinstance(arg,float):
1799          return math.sinh(arg)
1800       elif isinstance(arg,int):
1801          return math.sinh(arg)
1802       elif isinstance(arg,Symbol):
1803          return Sinh_Symbol(arg)
1804       else:
1805          raise TypeError,"sinh: Unknown argument type."
1806    
1807    class Sinh_Symbol(DependendSymbol):
1808       """
1809       L{Symbol} representing the result of the hyperbolic sine function
1810       """
1811       def __init__(self,arg):
1812          """
1813          initialization of sinh L{Symbol} with argument arg
1814          @param arg: argument of function
1815          @type arg: typically L{Symbol}.
1816          """
1817          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1818    
1819       def getMyCode(self,argstrs,format="escript"):
1820          """
1821          returns a program code that can be used to evaluate the symbol.
1822    
1823          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1824          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1825          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1826          @type format: C{str}
1827          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1828          @rtype: C{str}
1829          @raise: NotImplementedError: if the requested format is not available
1830          """
1831          if isinstance(argstrs,list):
1832              argstrs=argstrs[0]
1833          if format=="escript" or format=="str"  or format=="text":
1834             return "sinh(%s)"%argstrs
1835          else:
1836             raise NotImplementedError,"Sinh_Symbol does not provide program code for format %s."%format
1837    
1838       def substitute(self,argvals):
1839          """
1840          assigns new values to symbols in the definition of the symbol.
1841          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1842    
1843          @param argvals: new values assigned to symbols
1844          @type argvals: C{dict} with keywords of type L{Symbol}.
1845          @return: result of the substitution process. Operations are executed as much as possible.
1846          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1847          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1848          """
1849          if argvals.has_key(self):
1850             arg=argvals[self]
1851             if self.isAppropriateValue(arg):
1852                return arg
1853             else:
1854                raise TypeError,"%s: new value is not appropriate."%str(self)
1855          else:
1856             arg=self.getSubstitutedArguments(argvals)[0]
1857             return sinh(arg)
1858    
1859       def diff(self,arg):
1860          """
1861          differential of this object
1862    
1863          @param arg: the derivative is calculated with respect to arg
1864          @type arg: L{escript.Symbol}
1865          @return: derivative with respect to C{arg}
1866          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1867          """
1868          if arg==self:
1869             return identity(self.getShape())
1870          else:
1871             myarg=self.getArgument()[0]
1872             val=matchShape(cosh(myarg),self.getDifferentiatedArguments(arg)[0])
1873             return val[0]*val[1]
1874    
1875    def cosh(arg):
1876       """
1877       returns hyperbolic cosine of argument arg
1878    
1879       @param arg: argument
1880       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1881       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1882       @raises TypeError: if the type of the argument is not expected.
1883       """
1884       if isinstance(arg,numarray.NumArray):
1885          return numarray.cosh(arg)
1886       elif isinstance(arg,escript.Data):
1887          return arg._cosh()
1888       elif isinstance(arg,float):
1889          return math.cosh(arg)
1890       elif isinstance(arg,int):
1891          return math.cosh(arg)
1892       elif isinstance(arg,Symbol):
1893          return Cosh_Symbol(arg)
1894       else:
1895          raise TypeError,"cosh: Unknown argument type."
1896    
1897    class Cosh_Symbol(DependendSymbol):
1898       """
1899       L{Symbol} representing the result of the hyperbolic cosine function
1900       """
1901       def __init__(self,arg):
1902          """
1903          initialization of cosh L{Symbol} with argument arg
1904          @param arg: argument of function
1905          @type arg: typically L{Symbol}.
1906          """
1907          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1908    
1909       def getMyCode(self,argstrs,format="escript"):
1910          """
1911          returns a program code that can be used to evaluate the symbol.
1912    
1913          @param argstrs: gives for each argument a string representing the argument for the evaluation.
1914          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
1915          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
1916          @type format: C{str}
1917          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
1918          @rtype: C{str}
1919          @raise: NotImplementedError: if the requested format is not available
1920          """
1921          if isinstance(argstrs,list):
1922              argstrs=argstrs[0]
1923          if format=="escript" or format=="str"  or format=="text":
1924             return "cosh(%s)"%argstrs
1925          else:
1926             raise NotImplementedError,"Cosh_Symbol does not provide program code for format %s."%format
1927    
1928       def substitute(self,argvals):
1929          """
1930          assigns new values to symbols in the definition of the symbol.
1931          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
1932    
1933          @param argvals: new values assigned to symbols
1934          @type argvals: C{dict} with keywords of type L{Symbol}.
1935          @return: result of the substitution process. Operations are executed as much as possible.
1936          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
1937          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
1938          """
1939          if argvals.has_key(self):
1940             arg=argvals[self]
1941             if self.isAppropriateValue(arg):
1942                return arg
1943             else:
1944                raise TypeError,"%s: new value is not appropriate."%str(self)
1945          else:
1946             arg=self.getSubstitutedArguments(argvals)[0]
1947             return cosh(arg)
1948    
1949       def diff(self,arg):
1950          """
1951          differential of this object
1952    
1953          @param arg: the derivative is calculated with respect to arg
1954          @type arg: L{escript.Symbol}
1955          @return: derivative with respect to C{arg}
1956          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
1957          """
1958          if arg==self:
1959             return identity(self.getShape())
1960          else:
1961             myarg=self.getArgument()[0]
1962             val=matchShape(sinh(myarg),self.getDifferentiatedArguments(arg)[0])
1963             return val[0]*val[1]
1964    
1965    def tanh(arg):
1966       """
1967       returns hyperbolic tangent of argument arg
1968    
1969       @param arg: argument
1970       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1971       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1972       @raises TypeError: if the type of the argument is not expected.
1973       """
1974       if isinstance(arg,numarray.NumArray):
1975          return numarray.tanh(arg)
1976       elif isinstance(arg,escript.Data):
1977          return arg._tanh()
1978       elif isinstance(arg,float):
1979          return math.tanh(arg)
1980       elif isinstance(arg,int):
1981          return math.tanh(arg)
1982       elif isinstance(arg,Symbol):
1983          return Tanh_Symbol(arg)
1984       else:
1985          raise TypeError,"tanh: Unknown argument type."
1986    
1987    class Tanh_Symbol(DependendSymbol):
1988       """
1989       L{Symbol} representing the result of the hyperbolic tangent function
1990       """
1991       def __init__(self,arg):
1992          """
1993          initialization of tanh L{Symbol} with argument arg
1994          @param arg: argument of function
1995          @type arg: typically L{Symbol}.
1996          """
1997          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
1998    
1999       def getMyCode(self,argstrs,format="escript"):
2000          """
2001          returns a program code that can be used to evaluate the symbol.
2002    
2003          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2004          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2005          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2006          @type format: C{str}
2007          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2008          @rtype: C{str}
2009          @raise: NotImplementedError: if the requested format is not available
2010          """
2011          if isinstance(argstrs,list):
2012              argstrs=argstrs[0]
2013          if format=="escript" or format=="str"  or format=="text":
2014             return "tanh(%s)"%argstrs
2015          else:
2016             raise NotImplementedError,"Tanh_Symbol does not provide program code for format %s."%format
2017    
2018       def substitute(self,argvals):
2019          """
2020          assigns new values to symbols in the definition of the symbol.
2021          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2022    
2023          @param argvals: new values assigned to symbols
2024          @type argvals: C{dict} with keywords of type L{Symbol}.
2025          @return: result of the substitution process. Operations are executed as much as possible.
2026          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2027          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2028          """
2029          if argvals.has_key(self):
2030             arg=argvals[self]
2031             if self.isAppropriateValue(arg):
2032                return arg
2033             else:
2034                raise TypeError,"%s: new value is not appropriate."%str(self)
2035          else:
2036             arg=self.getSubstitutedArguments(argvals)[0]
2037             return tanh(arg)
2038    
2039       def diff(self,arg):
2040          """
2041          differential of this object
2042    
2043          @param arg: the derivative is calculated with respect to arg
2044          @type arg: L{escript.Symbol}
2045          @return: derivative with respect to C{arg}
2046          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2047          """
2048          if arg==self:
2049             return identity(self.getShape())
2050          else:
2051             myarg=self.getArgument()[0]
2052             val=matchShape(1./cosh(myarg)**2,self.getDifferentiatedArguments(arg)[0])
2053             return val[0]*val[1]
2054    
2055    def asinh(arg):
2056       """
2057       returns inverse hyperbolic sine of argument arg
2058    
2059       @param arg: argument
2060       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2061       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2062       @raises TypeError: if the type of the argument is not expected.
2063       """
2064       if isinstance(arg,numarray.NumArray):
2065          return numarray.arcsinh(arg)
2066       elif isinstance(arg,escript.Data):
2067          return arg._asinh()
2068       elif isinstance(arg,float):
2069          return numarray.arcsinh(arg)
2070       elif isinstance(arg,int):
2071          return numarray.arcsinh(float(arg))
2072       elif isinstance(arg,Symbol):
2073          return Asinh_Symbol(arg)
2074       else:
2075          raise TypeError,"asinh: Unknown argument type."
2076    
2077    class Asinh_Symbol(DependendSymbol):
2078       """
2079       L{Symbol} representing the result of the inverse hyperbolic sine function
2080       """
2081       def __init__(self,arg):
2082          """
2083          initialization of asinh L{Symbol} with argument arg
2084          @param arg: argument of function
2085          @type arg: typically L{Symbol}.
2086          """
2087          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2088    
2089       def getMyCode(self,argstrs,format="escript"):
2090          """
2091          returns a program code that can be used to evaluate the symbol.
2092    
2093          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2094          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2095          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2096          @type format: C{str}
2097          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2098          @rtype: C{str}
2099          @raise: NotImplementedError: if the requested format is not available
2100          """
2101          if isinstance(argstrs,list):
2102              argstrs=argstrs[0]
2103          if format=="escript" or format=="str"  or format=="text":
2104             return "asinh(%s)"%argstrs
2105          else:
2106             raise NotImplementedError,"Asinh_Symbol does not provide program code for format %s."%format
2107    
2108       def substitute(self,argvals):
2109          """
2110          assigns new values to symbols in the definition of the symbol.
2111          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2112    
2113          @param argvals: new values assigned to symbols
2114          @type argvals: C{dict} with keywords of type L{Symbol}.
2115          @return: result of the substitution process. Operations are executed as much as possible.
2116          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2117          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2118          """
2119          if argvals.has_key(self):
2120             arg=argvals[self]
2121             if self.isAppropriateValue(arg):
2122                return arg
2123             else:
2124                raise TypeError,"%s: new value is not appropriate."%str(self)
2125          else:
2126             arg=self.getSubstitutedArguments(argvals)[0]
2127             return asinh(arg)
2128    
2129       def diff(self,arg):
2130          """
2131          differential of this object
2132    
2133          @param arg: the derivative is calculated with respect to arg
2134          @type arg: L{escript.Symbol}
2135          @return: derivative with respect to C{arg}
2136          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2137          """
2138          if arg==self:
2139             return identity(self.getShape())
2140          else:
2141             myarg=self.getArgument()[0]
2142             val=matchShape(1./sqrt(myarg**2+1),self.getDifferentiatedArguments(arg)[0])
2143             return val[0]*val[1]
2144    
2145    def acosh(arg):
2146       """
2147       returns inverse hyperolic cosine of argument arg
2148    
2149       @param arg: argument
2150       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2151       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2152       @raises TypeError: if the type of the argument is not expected.
2153       """
2154       if isinstance(arg,numarray.NumArray):
2155          return numarray.arccosh(arg)
2156       elif isinstance(arg,escript.Data):
2157          return arg._acosh()
2158       elif isinstance(arg,float):
2159          return numarray.arccosh(arg)
2160       elif isinstance(arg,int):
2161          return numarray.arccosh(float(arg))
2162       elif isinstance(arg,Symbol):
2163          return Acosh_Symbol(arg)
2164       else:
2165          raise TypeError,"acosh: Unknown argument type."
2166    
2167    class Acosh_Symbol(DependendSymbol):
2168       """
2169       L{Symbol} representing the result of the inverse hyperolic cosine function
2170       """
2171       def __init__(self,arg):
2172          """
2173          initialization of acosh L{Symbol} with argument arg
2174          @param arg: argument of function
2175          @type arg: typically L{Symbol}.
2176          """
2177          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2178    
2179       def getMyCode(self,argstrs,format="escript"):
2180          """
2181          returns a program code that can be used to evaluate the symbol.
2182    
2183          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2184          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2185          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2186          @type format: C{str}
2187          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2188          @rtype: C{str}
2189          @raise: NotImplementedError: if the requested format is not available
2190          """
2191          if isinstance(argstrs,list):
2192              argstrs=argstrs[0]
2193          if format=="escript" or format=="str"  or format=="text":
2194             return "acosh(%s)"%argstrs
2195          else:
2196             raise NotImplementedError,"Acosh_Symbol does not provide program code for format %s."%format
2197    
2198       def substitute(self,argvals):
2199          """
2200          assigns new values to symbols in the definition of the symbol.
2201          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2202    
2203          @param argvals: new values assigned to symbols
2204          @type argvals: C{dict} with keywords of type L{Symbol}.
2205          @return: result of the substitution process. Operations are executed as much as possible.
2206          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2207          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2208          """
2209          if argvals.has_key(self):
2210             arg=argvals[self]
2211             if self.isAppropriateValue(arg):
2212                return arg
2213             else:
2214                raise TypeError,"%s: new value is not appropriate."%str(self)
2215          else:
2216             arg=self.getSubstitutedArguments(argvals)[0]
2217             return acosh(arg)
2218    
2219       def diff(self,arg):
2220          """
2221          differential of this object
2222    
2223          @param arg: the derivative is calculated with respect to arg
2224          @type arg: L{escript.Symbol}
2225          @return: derivative with respect to C{arg}
2226          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2227          """
2228          if arg==self:
2229             return identity(self.getShape())
2230          else:
2231             myarg=self.getArgument()[0]
2232             val=matchShape(1./sqrt(myarg**2-1),self.getDifferentiatedArguments(arg)[0])
2233             return val[0]*val[1]
2234    
2235    def atanh(arg):
2236       """
2237       returns inverse hyperbolic tangent of argument arg
2238    
2239       @param arg: argument
2240       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2241       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2242       @raises TypeError: if the type of the argument is not expected.
2243       """
2244       if isinstance(arg,numarray.NumArray):
2245          return numarray.arctanh(arg)
2246       elif isinstance(arg,escript.Data):
2247          return arg._atanh()
2248       elif isinstance(arg,float):
2249          return numarray.arctanh(arg)
2250       elif isinstance(arg,int):
2251          return numarray.arctanh(float(arg))
2252       elif isinstance(arg,Symbol):
2253          return Atanh_Symbol(arg)
2254       else:
2255          raise TypeError,"atanh: Unknown argument type."
2256    
2257    class Atanh_Symbol(DependendSymbol):
2258       """
2259       L{Symbol} representing the result of the inverse hyperbolic tangent function
2260       """
2261       def __init__(self,arg):
2262          """
2263          initialization of atanh L{Symbol} with argument arg
2264          @param arg: argument of function
2265          @type arg: typically L{Symbol}.
2266          """
2267          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2268    
2269       def getMyCode(self,argstrs,format="escript"):
2270          """
2271          returns a program code that can be used to evaluate the symbol.
2272    
2273          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2274          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2275          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2276          @type format: C{str}
2277          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2278          @rtype: C{str}
2279          @raise: NotImplementedError: if the requested format is not available
2280          """
2281          if isinstance(argstrs,list):
2282              argstrs=argstrs[0]
2283          if format=="escript" or format=="str"  or format=="text":
2284             return "atanh(%s)"%argstrs
2285          else:
2286             raise NotImplementedError,"Atanh_Symbol does not provide program code for format %s."%format
2287    
2288       def substitute(self,argvals):
2289          """
2290          assigns new values to symbols in the definition of the symbol.
2291          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2292    
2293          @param argvals: new values assigned to symbols
2294          @type argvals: C{dict} with keywords of type L{Symbol}.
2295          @return: result of the substitution process. Operations are executed as much as possible.
2296          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2297          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2298          """
2299          if argvals.has_key(self):
2300             arg=argvals[self]
2301             if self.isAppropriateValue(arg):
2302                return arg
2303             else:
2304                raise TypeError,"%s: new value is not appropriate."%str(self)
2305          else:
2306             arg=self.getSubstitutedArguments(argvals)[0]
2307             return atanh(arg)
2308    
2309       def diff(self,arg):
2310          """
2311          differential of this object
2312    
2313          @param arg: the derivative is calculated with respect to arg
2314          @type arg: L{escript.Symbol}
2315          @return: derivative with respect to C{arg}
2316          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2317          """
2318          if arg==self:
2319             return identity(self.getShape())
2320          else:
2321             myarg=self.getArgument()[0]
2322             val=matchShape(1./(1.-myarg**2),self.getDifferentiatedArguments(arg)[0])
2323             return val[0]*val[1]
2324    
2325    def exp(arg):
2326       """
2327       returns exponential of argument arg
2328    
2329       @param arg: argument
2330       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2331       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2332       @raises TypeError: if the type of the argument is not expected.
2333       """
2334       if isinstance(arg,numarray.NumArray):
2335          return numarray.exp(arg)
2336       elif isinstance(arg,escript.Data):
2337          return arg._exp()
2338       elif isinstance(arg,float):
2339          return math.exp(arg)
2340       elif isinstance(arg,int):
2341          return math.exp(arg)
2342       elif isinstance(arg,Symbol):
2343          return Exp_Symbol(arg)
2344       else:
2345          raise TypeError,"exp: Unknown argument type."
2346    
2347    class Exp_Symbol(DependendSymbol):
2348       """
2349       L{Symbol} representing the result of the exponential function
2350       """
2351       def __init__(self,arg):
2352          """
2353          initialization of exp L{Symbol} with argument arg
2354          @param arg: argument of function
2355          @type arg: typically L{Symbol}.
2356          """
2357          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2358    
2359       def getMyCode(self,argstrs,format="escript"):
2360          """
2361          returns a program code that can be used to evaluate the symbol.
2362    
2363          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2364          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2365          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2366          @type format: C{str}
2367          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2368          @rtype: C{str}
2369          @raise: NotImplementedError: if the requested format is not available
2370          """
2371          if isinstance(argstrs,list):
2372              argstrs=argstrs[0]
2373          if format=="escript" or format=="str"  or format=="text":
2374             return "exp(%s)"%argstrs
2375          else:
2376             raise NotImplementedError,"Exp_Symbol does not provide program code for format %s."%format
2377    
2378       def substitute(self,argvals):
2379          """
2380          assigns new values to symbols in the definition of the symbol.
2381          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2382    
2383          @param argvals: new values assigned to symbols
2384          @type argvals: C{dict} with keywords of type L{Symbol}.
2385          @return: result of the substitution process. Operations are executed as much as possible.
2386          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2387          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2388          """
2389          if argvals.has_key(self):
2390             arg=argvals[self]
2391             if self.isAppropriateValue(arg):
2392                return arg
2393             else:
2394                raise TypeError,"%s: new value is not appropriate."%str(self)
2395          else:
2396             arg=self.getSubstitutedArguments(argvals)[0]
2397             return exp(arg)
2398    
2399       def diff(self,arg):
2400          """
2401          differential of this object
2402    
2403          @param arg: the derivative is calculated with respect to arg
2404          @type arg: L{escript.Symbol}
2405          @return: derivative with respect to C{arg}
2406          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2407          """
2408          if arg==self:
2409             return identity(self.getShape())
2410          else:
2411             myarg=self.getArgument()[0]
2412             val=matchShape(self,self.getDifferentiatedArguments(arg)[0])
2413             return val[0]*val[1]
2414    
2415    def sqrt(arg):
2416       """
2417       returns square root of argument arg
2418    
2419       @param arg: argument
2420       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2421       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2422       @raises TypeError: if the type of the argument is not expected.
2423       """
2424       if isinstance(arg,numarray.NumArray):
2425          return numarray.sqrt(arg)
2426       elif isinstance(arg,escript.Data):
2427          return arg._sqrt()
2428       elif isinstance(arg,float):
2429          return math.sqrt(arg)
2430       elif isinstance(arg,int):
2431          return math.sqrt(arg)
2432       elif isinstance(arg,Symbol):
2433          return Sqrt_Symbol(arg)
2434       else:
2435          raise TypeError,"sqrt: Unknown argument type."
2436    
2437    class Sqrt_Symbol(DependendSymbol):
2438       """
2439       L{Symbol} representing the result of the square root function
2440       """
2441       def __init__(self,arg):
2442          """
2443          initialization of sqrt L{Symbol} with argument arg
2444          @param arg: argument of function
2445          @type arg: typically L{Symbol}.
2446          """
2447          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2448    
2449       def getMyCode(self,argstrs,format="escript"):
2450          """
2451          returns a program code that can be used to evaluate the symbol.
2452    
2453          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2454          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2455          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2456          @type format: C{str}
2457          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2458          @rtype: C{str}
2459          @raise: NotImplementedError: if the requested format is not available
2460          """
2461          if isinstance(argstrs,list):
2462              argstrs=argstrs[0]
2463          if format=="escript" or format=="str"  or format=="text":
2464             return "sqrt(%s)"%argstrs
2465          else:
2466             raise NotImplementedError,"Sqrt_Symbol does not provide program code for format %s."%format
2467    
2468       def substitute(self,argvals):
2469          """
2470          assigns new values to symbols in the definition of the symbol.
2471          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2472    
2473          @param argvals: new values assigned to symbols
2474          @type argvals: C{dict} with keywords of type L{Symbol}.
2475          @return: result of the substitution process. Operations are executed as much as possible.
2476          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2477          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2478          """
2479          if argvals.has_key(self):
2480             arg=argvals[self]
2481             if self.isAppropriateValue(arg):
2482                return arg
2483             else:
2484                raise TypeError,"%s: new value is not appropriate."%str(self)
2485          else:
2486             arg=self.getSubstitutedArguments(argvals)[0]
2487             return sqrt(arg)
2488    
2489       def diff(self,arg):
2490          """
2491          differential of this object
2492    
2493          @param arg: the derivative is calculated with respect to arg
2494          @type arg: L{escript.Symbol}
2495          @return: derivative with respect to C{arg}
2496          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2497          """
2498          if arg==self:
2499             return identity(self.getShape())
2500          else:
2501             myarg=self.getArgument()[0]
2502             val=matchShape(0.5/self,self.getDifferentiatedArguments(arg)[0])
2503             return val[0]*val[1]
2504    
2505    def log(arg):
2506       """
2507       returns natural logarithm of argument arg
2508    
2509       @param arg: argument
2510       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2511       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2512       @raises TypeError: if the type of the argument is not expected.
2513       """
2514       if isinstance(arg,numarray.NumArray):
2515          return numarray.log(arg)
2516       elif isinstance(arg,escript.Data):
2517          return arg._log()
2518       elif isinstance(arg,float):
2519          return math.log(arg)
2520       elif isinstance(arg,int):
2521          return math.log(arg)
2522       elif isinstance(arg,Symbol):
2523          return Log_Symbol(arg)
2524       else:
2525          raise TypeError,"log: Unknown argument type."
2526    
2527    class Log_Symbol(DependendSymbol):
2528       """
2529       L{Symbol} representing the result of the natural logarithm function
2530       """
2531       def __init__(self,arg):
2532          """
2533          initialization of log L{Symbol} with argument arg
2534          @param arg: argument of function
2535          @type arg: typically L{Symbol}.
2536          """
2537          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2538    
2539       def getMyCode(self,argstrs,format="escript"):
2540          """
2541          returns a program code that can be used to evaluate the symbol.
2542    
2543          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2544          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2545          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2546          @type format: C{str}
2547          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2548          @rtype: C{str}
2549          @raise: NotImplementedError: if the requested format is not available
2550          """
2551          if isinstance(argstrs,list):
2552              argstrs=argstrs[0]
2553          if format=="escript" or format=="str"  or format=="text":
2554             return "log(%s)"%argstrs
2555          else:
2556             raise NotImplementedError,"Log_Symbol does not provide program code for format %s."%format
2557    
2558       def substitute(self,argvals):
2559          """
2560          assigns new values to symbols in the definition of the symbol.
2561          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2562    
2563          @param argvals: new values assigned to symbols
2564          @type argvals: C{dict} with keywords of type L{Symbol}.
2565          @return: result of the substitution process. Operations are executed as much as possible.
2566          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2567          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2568          """
2569          if argvals.has_key(self):
2570             arg=argvals[self]
2571             if self.isAppropriateValue(arg):
2572                return arg
2573             else:
2574                raise TypeError,"%s: new value is not appropriate."%str(self)
2575          else:
2576             arg=self.getSubstitutedArguments(argvals)[0]
2577             return log(arg)
2578    
2579       def diff(self,arg):
2580          """
2581          differential of this object
2582    
2583          @param arg: the derivative is calculated with respect to arg
2584          @type arg: L{escript.Symbol}
2585          @return: derivative with respect to C{arg}
2586          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2587          """
2588          if arg==self:
2589             return identity(self.getShape())
2590          else:
2591             myarg=self.getArgument()[0]
2592             val=matchShape(1./arg,self.getDifferentiatedArguments(arg)[0])
2593             return val[0]*val[1]
2594    
2595    def sign(arg):
2596       """
2597       returns sign of argument arg
2598    
2599       @param arg: argument
2600       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2601       @rtype:C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2602       @raises TypeError: if the type of the argument is not expected.
2603       """
2604       if isinstance(arg,numarray.NumArray):
2605          return wherePositive(arg)-whereNegative(arg)
2606       elif isinstance(arg,escript.Data):
2607          return arg._sign()
2608       elif isinstance(arg,float):
2609          if arg>0:
2610            return 1.
2611          elif arg<0:
2612            return -1.
2613          else:
2614            return 0.
2615       elif isinstance(arg,int):
2616          if float(arg)>0:
2617            return 1.
2618          elif float(arg)<0:
2619            return -1.
2620          else:
2621            return 0.
2622       elif isinstance(arg,Symbol):
2623          return wherePositive(arg)-whereNegative(arg)
2624       else:
2625          raise TypeError,"sign: Unknown argument type."
2626    
2627    class Abs_Symbol(DependendSymbol):
2628       """
2629       L{Symbol} representing the result of the absolute value function
2630       """
2631       def __init__(self,arg):
2632          """
2633          initialization of abs L{Symbol} with argument arg
2634          @param arg: argument of function
2635          @type arg: typically L{Symbol}.
2636          """
2637          DependendSymbol.__init__(self,args=[arg],shape=arg.getShape(),dim=arg.getDim())
2638    
2639       def getMyCode(self,argstrs,format="escript"):
2640          """
2641          returns a program code that can be used to evaluate the symbol.
2642    
2643          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2644          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2645          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2646          @type format: C{str}
2647          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2648          @rtype: C{str}
2649          @raise: NotImplementedError: if the requested format is not available
2650          """
2651          if isinstance(argstrs,list):
2652              argstrs=argstrs[0]
2653          if format=="escript" or format=="str"  or format=="text":
2654             return "abs(%s)"%argstrs
2655          else:
2656             raise NotImplementedError,"Abs_Symbol does not provide program code for format %s."%format
2657    
2658       def substitute(self,argvals):
2659          """
2660          assigns new values to symbols in the definition of the symbol.
2661          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2662    
2663          @param argvals: new values assigned to symbols
2664          @type argvals: C{dict} with keywords of type L{Symbol}.
2665          @return: result of the substitution process. Operations are executed as much as possible.
2666          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2667          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2668          """
2669          if argvals.has_key(self):
2670             arg=argvals[self]
2671             if self.isAppropriateValue(arg):
2672                return arg
2673             else:
2674                raise TypeError,"%s: new value is not appropriate."%str(self)
2675          else:
2676             arg=self.getSubstitutedArguments(argvals)[0]
2677             return abs(arg)
2678    
2679       def diff(self,arg):
2680          """
2681          differential of this object
2682    
2683          @param arg: the derivative is calculated with respect to arg
2684          @type arg: L{escript.Symbol}
2685          @return: derivative with respect to C{arg}
2686          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
2687          """
2688          if arg==self:
2689             return identity(self.getShape())
2690          else:
2691             myarg=self.getArgument()[0]
2692             val=matchShape(sign(myarg),self.getDifferentiatedArguments(arg)[0])
2693             return val[0]*val[1]
2694    
2695    def minval(arg):
2696       """
2697       returns minimum value over all components of arg at each data point
2698    
2699       @param arg: argument
2700       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2701       @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2702       @raises TypeError: if the type of the argument is not expected.
2703       """
2704       if isinstance(arg,numarray.NumArray):
2705          if arg.rank==0:
2706             return float(arg)
2707          else:
2708             return arg.min()
2709       elif isinstance(arg,escript.Data):
2710          return arg._minval()
2711       elif isinstance(arg,float):
2712          return arg
2713       elif isinstance(arg,int):
2714          return float(arg)
2715       elif isinstance(arg,Symbol):
2716          return Minval_Symbol(arg)
2717       else:
2718          raise TypeError,"minval: Unknown argument type."
2719    
2720    class Minval_Symbol(DependendSymbol):
2721       """
2722       L{Symbol} representing the result of the minimum value function
2723       """
2724       def __init__(self,arg):
2725          """
2726          initialization of minimum value L{Symbol} with argument arg
2727          @param arg: argument of function
2728          @type arg: typically L{Symbol}.
2729          """
2730          DependendSymbol.__init__(self,args=[arg],shape=(),dim=arg.getDim())
2731    
2732       def getMyCode(self,argstrs,format="escript"):
2733          """
2734          returns a program code that can be used to evaluate the symbol.
2735    
2736          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2737          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2738          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2739          @type format: C{str}
2740          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2741          @rtype: C{str}
2742          @raise: NotImplementedError: if the requested format is not available
2743          """
2744          if isinstance(argstrs,list):
2745              argstrs=argstrs[0]
2746          if format=="escript" or format=="str"  or format=="text":
2747             return "minval(%s)"%argstrs
2748          else:
2749             raise NotImplementedError,"Minval_Symbol does not provide program code for format %s."%format
2750    
2751       def substitute(self,argvals):
2752          """
2753          assigns new values to symbols in the definition of the symbol.
2754          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2755    
2756          @param argvals: new values assigned to symbols
2757          @type argvals: C{dict} with keywords of type L{Symbol}.
2758          @return: result of the substitution process. Operations are executed as much as possible.
2759          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2760          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2761          """
2762          if argvals.has_key(self):
2763             arg=argvals[self]
2764             if self.isAppropriateValue(arg):
2765                return arg
2766             else:
2767                raise TypeError,"%s: new value is not appropriate."%str(self)
2768          else:
2769             arg=self.getSubstitutedArguments(argvals)[0]
2770             return minval(arg)
2771    
2772    def maxval(arg):
2773       """
2774       returns maximum value over all components of arg at each data point
2775    
2776       @param arg: argument
2777       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2778       @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2779       @raises TypeError: if the type of the argument is not expected.
2780       """
2781       if isinstance(arg,numarray.NumArray):
2782          if arg.rank==0:
2783             return float(arg)
2784          else:
2785             return arg.max()
2786       elif isinstance(arg,escript.Data):
2787          return arg._maxval()
2788       elif isinstance(arg,float):
2789          return arg
2790       elif isinstance(arg,int):
2791          return float(arg)
2792       elif isinstance(arg,Symbol):
2793          return Maxval_Symbol(arg)
2794       else:
2795          raise TypeError,"maxval: Unknown argument type."
2796    
2797    class Maxval_Symbol(DependendSymbol):
2798       """
2799       L{Symbol} representing the result of the maximum value function
2800       """
2801       def __init__(self,arg):
2802          """
2803          initialization of maximum value L{Symbol} with argument arg
2804          @param arg: argument of function
2805          @type arg: typically L{Symbol}.
2806          """
2807          DependendSymbol.__init__(self,args=[arg],shape=(),dim=arg.getDim())
2808    
2809       def getMyCode(self,argstrs,format="escript"):
2810          """
2811          returns a program code that can be used to evaluate the symbol.
2812    
2813          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2814          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2815          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2816          @type format: C{str}
2817          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2818          @rtype: C{str}
2819          @raise: NotImplementedError: if the requested format is not available
2820          """
2821          if isinstance(argstrs,list):
2822              argstrs=argstrs[0]
2823          if format=="escript" or format=="str"  or format=="text":
2824             return "maxval(%s)"%argstrs
2825          else:
2826             raise NotImplementedError,"Maxval_Symbol does not provide program code for format %s."%format
2827    
2828       def substitute(self,argvals):
2829          """
2830          assigns new values to symbols in the definition of the symbol.
2831          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2832    
2833          @param argvals: new values assigned to symbols
2834          @type argvals: C{dict} with keywords of type L{Symbol}.
2835          @return: result of the substitution process. Operations are executed as much as possible.
2836          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2837          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2838          """
2839          if argvals.has_key(self):
2840             arg=argvals[self]
2841             if self.isAppropriateValue(arg):
2842                return arg
2843             else:
2844                raise TypeError,"%s: new value is not appropriate."%str(self)
2845          else:
2846             arg=self.getSubstitutedArguments(argvals)[0]
2847             return maxval(arg)
2848    
2849  def length(arg):  def length(arg):
2850       """
2851       returns length/Euclidean norm of argument arg at each data point
2852    
2853       @param arg: argument
2854       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2855       @rtype:C{float}, L{escript.Data}, L{Symbol} depending on the type of arg.
2856       """
2857       return sqrt(inner(arg,arg))
2858    
2859    def trace(arg,axis_offset=0):
2860       """
2861       returns the trace of arg which the sum of arg[k,k] over k.
2862    
2863       @param arg: argument
2864       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2865       @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2866                      axis_offset and axis_offset+1 must be equal.
2867       @type axis_offset: C{int}
2868       @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2869       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
2870       """
2871       if isinstance(arg,numarray.NumArray):
2872          sh=arg.shape
2873          if len(sh)<2:
2874            raise ValueError,"trace: rank of argument must be greater than 1"
2875          if axis_offset<0 or axis_offset>len(sh)-2:
2876            raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2
2877          s1=1
2878          for i in range(axis_offset): s1*=sh[i]
2879          s2=1
2880          for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2881          if not sh[axis_offset] == sh[axis_offset+1]:
2882            raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2883          arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2884          out=numarray.zeros([s1,s2],numarray.Float)
2885          for i1 in range(s1):
2886            for i2 in range(s2):
2887                for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
2888          out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2889          return out
2890       elif isinstance(arg,escript.Data):
2891          return escript_trace(arg,axis_offset)
2892       elif isinstance(arg,float):
2893          raise TypeError,"trace: illegal argument type float."
2894       elif isinstance(arg,int):
2895          raise TypeError,"trace: illegal argument type int."
2896       elif isinstance(arg,Symbol):
2897          return Trace_Symbol(arg,axis_offset)
2898       else:
2899          raise TypeError,"trace: Unknown argument type."
2900    
2901    def escript_trace(arg,axis_offset): # this should be escript._trace
2902          "arg si a Data objects!!!"
2903          if arg.getRank()<2:
2904            raise ValueError,"escript_trace: rank of argument must be greater than 1"
2905          if axis_offset<0 or axis_offset>arg.getRank()-2:
2906            raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2
2907          s=list(arg.getShape())        
2908          if not s[axis_offset] == s[axis_offset+1]:
2909            raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2910          out=escript.Data(0.,tuple(s[0:axis_offset]+s[axis_offset+2:]),arg.getFunctionSpace())
2911          if arg.getRank()==2:
2912             for i0 in range(s[0]):
2913                out+=arg[i0,i0]
2914          elif arg.getRank()==3:
2915             if axis_offset==0:
2916                for i0 in range(s[0]):
2917                      for i2 in range(s[2]):
2918                             out[i2]+=arg[i0,i0,i2]
2919             elif axis_offset==1:
2920                for i0 in range(s[0]):
2921                   for i1 in range(s[1]):
2922                             out[i0]+=arg[i0,i1,i1]
2923          elif arg.getRank()==4:
2924             if axis_offset==0:
2925                for i0 in range(s[0]):
2926                      for i2 in range(s[2]):
2927                         for i3 in range(s[3]):
2928                             out[i2,i3]+=arg[i0,i0,i2,i3]
2929             elif axis_offset==1:
2930                for i0 in range(s[0]):
2931                   for i1 in range(s[1]):
2932                         for i3 in range(s[3]):
2933                             out[i0,i3]+=arg[i0,i1,i1,i3]
2934             elif axis_offset==2:
2935                for i0 in range(s[0]):
2936                   for i1 in range(s[1]):
2937                      for i2 in range(s[2]):
2938                             out[i0,i1]+=arg[i0,i1,i2,i2]
2939          return out
2940    class Trace_Symbol(DependendSymbol):
2941       """
2942       L{Symbol} representing the result of the trace function
2943       """
2944       def __init__(self,arg,axis_offset=0):
2945          """
2946          initialization of trace L{Symbol} with argument arg
2947          @param arg: argument of function
2948          @type arg: L{Symbol}.
2949          @param axis_offset: axis_offset to components to sum over. C{axis_offset} must be non-negative and less than the rank of arg +1. The dimensions on component
2950                      axis_offset and axis_offset+1 must be equal.
2951          @type axis_offset: C{int}
2952          """
2953          if arg.getRank()<2:
2954            raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"
2955          if axis_offset<0 or axis_offset>arg.getRank()-2:
2956            raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2
2957          s=list(arg.getShape())        
2958          if not s[axis_offset] == s[axis_offset+1]:
2959            raise ValueError,"Trace_Symbol: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2960          super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
2961    
2962       def getMyCode(self,argstrs,format="escript"):
2963          """
2964          returns a program code that can be used to evaluate the symbol.
2965    
2966          @param argstrs: gives for each argument a string representing the argument for the evaluation.
2967          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
2968          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
2969          @type format: C{str}
2970          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
2971          @rtype: C{str}
2972          @raise: NotImplementedError: if the requested format is not available
2973          """
2974          if format=="escript" or format=="str"  or format=="text":
2975             return "trace(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
2976          else:
2977             raise NotImplementedError,"Trace_Symbol does not provide program code for format %s."%format
2978    
2979       def substitute(self,argvals):
2980          """
2981          assigns new values to symbols in the definition of the symbol.
2982          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
2983    
2984          @param argvals: new values assigned to symbols
2985          @type argvals: C{dict} with keywords of type L{Symbol}.
2986          @return: result of the substitution process. Operations are executed as much as possible.
2987          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
2988          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
2989          """
2990          if argvals.has_key(self):
2991             arg=argvals[self]
2992             if self.isAppropriateValue(arg):
2993                return arg
2994             else:
2995                raise TypeError,"%s: new value is not appropriate."%str(self)
2996          else:
2997             arg=self.getSubstitutedArguments(argvals)
2998             return trace(arg[0],axis_offset=arg[1])
2999    
3000       def diff(self,arg):
3001          """
3002          differential of this object
3003    
3004          @param arg: the derivative is calculated with respect to arg
3005          @type arg: L{escript.Symbol}
3006          @return: derivative with respect to C{arg}
3007          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3008          """
3009          if arg==self:
3010             return identity(self.getShape())
3011          else:
3012             return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3013    
3014    def transpose(arg,axis_offset=None):
3015       """
3016       returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
3017    
3018       @param arg: argument
3019       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3020       @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3021                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3022       @type axis_offset: C{int}
3023       @return: transpose of arg
3024       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3025       """
3026       if isinstance(arg,numarray.NumArray):
3027          if axis_offset==None: axis_offset=int(arg.rank/2)
3028          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3029       elif isinstance(arg,escript.Data):
3030          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3031          return escript_transpose(arg,axis_offset)
3032       elif isinstance(arg,float):
3033          if not ( axis_offset==0 or axis_offset==None):
3034            raise ValueError,"transpose: axis_offset must be 0 for float argument"
3035          return arg
3036       elif isinstance(arg,int):
3037          if not ( axis_offset==0 or axis_offset==None):
3038            raise ValueError,"transpose: axis_offset must be 0 for int argument"
3039          return float(arg)
3040       elif isinstance(arg,Symbol):
3041          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3042          return Transpose_Symbol(arg,axis_offset)
3043       else:
3044          raise TypeError,"transpose: Unknown argument type."
3045    
3046    def escript_transpose(arg,axis_offset): # this should be escript._transpose
3047          "arg si a Data objects!!!"
3048          r=arg.getRank()
3049          if axis_offset<0 or axis_offset>r:
3050            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3051          s=arg.getShape()
3052          s_out=s[axis_offset:]+s[:axis_offset]
3053          out=escript.Data(0.,s_out,arg.getFunctionSpace())
3054          if r==4:
3055             if axis_offset==1:
3056                for i0 in range(s_out[0]):
3057                   for i1 in range(s_out[1]):
3058                      for i2 in range(s_out[2]):
3059                         for i3 in range(s_out[3]):
3060                             out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]
3061             elif axis_offset==2:
3062                for i0 in range(s_out[0]):
3063                   for i1 in range(s_out[1]):
3064                      for i2 in range(s_out[2]):
3065                         for i3 in range(s_out[3]):
3066                             out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]
3067             elif axis_offset==3:
3068                for i0 in range(s_out[0]):
3069                   for i1 in range(s_out[1]):
3070                      for i2 in range(s_out[2]):
3071                         for i3 in range(s_out[3]):
3072                             out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]
3073             else:
3074                for i0 in range(s_out[0]):
3075                   for i1 in range(s_out[1]):
3076                      for i2 in range(s_out[2]):
3077                         for i3 in range(s_out[3]):
3078                             out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]
3079          elif r==3:
3080             if axis_offset==1:
3081                for i0 in range(s_out[0]):
3082                   for i1 in range(s_out[1]):
3083                      for i2 in range(s_out[2]):
3084                             out[i0,i1,i2]=arg[i2,i0,i1]
3085             elif axis_offset==2:
3086                for i0 in range(s_out[0]):
3087                   for i1 in range(s_out[1]):
3088                      for i2 in range(s_out[2]):
3089                             out[i0,i1,i2]=arg[i1,i2,i0]
3090             else:
3091                for i0 in range(s_out[0]):
3092                   for i1 in range(s_out[1]):
3093                      for i2 in range(s_out[2]):
3094                             out[i0,i1,i2]=arg[i0,i1,i2]
3095          elif r==2:
3096             if axis_offset==1:
3097                for i0 in range(s_out[0]):
3098                   for i1 in range(s_out[1]):
3099                             out[i0,i1]=arg[i1,i0]
3100             else:
3101                for i0 in range(s_out[0]):
3102                   for i1 in range(s_out[1]):
3103                             out[i0,i1]=arg[i0,i1]
3104          elif r==1:
3105              for i0 in range(s_out[0]):
3106                   out[i0]=arg[i0]
3107          elif r==0:
3108                 out=arg+0.
3109          return out
3110    class Transpose_Symbol(DependendSymbol):
3111       """
3112       L{Symbol} representing the result of the transpose function
3113       """
3114       def __init__(self,arg,axis_offset=None):
3115          """
3116          initialization of transpose L{Symbol} with argument arg
3117    
3118          @param arg: argument of function
3119          @type arg: L{Symbol}.
3120           @param axis_offset: the first axis_offset components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3121                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3122          @type axis_offset: C{int}
3123          """
3124          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3125          if axis_offset<0 or axis_offset>arg.getRank():
3126            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3127          s=arg.getShape()
3128          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3129    
3130       def getMyCode(self,argstrs,format="escript"):
3131          """
3132          returns a program code that can be used to evaluate the symbol.
3133    
3134          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3135          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3136          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3137          @type format: C{str}
3138          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3139          @rtype: C{str}
3140          @raise: NotImplementedError: if the requested format is not available
3141          """
3142          if format=="escript" or format=="str"  or format=="text":
3143             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3144          else:
3145             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3146    
3147       def substitute(self,argvals):
3148          """
3149          assigns new values to symbols in the definition of the symbol.
3150          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3151    
3152          @param argvals: new values assigned to symbols
3153          @type argvals: C{dict} with keywords of type L{Symbol}.
3154          @return: result of the substitution process. Operations are executed as much as possible.
3155          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3156          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3157          """
3158          if argvals.has_key(self):
3159             arg=argvals[self]
3160             if self.isAppropriateValue(arg):
3161                return arg
3162             else:
3163                raise TypeError,"%s: new value is not appropriate."%str(self)
3164          else:
3165             arg=self.getSubstitutedArguments(argvals)
3166             return transpose(arg[0],axis_offset=arg[1])
3167    
3168       def diff(self,arg):
3169          """
3170          differential of this object
3171    
3172          @param arg: the derivative is calculated with respect to arg
3173          @type arg: L{escript.Symbol}
3174          @return: derivative with respect to C{arg}
3175          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3176          """
3177          if arg==self:
3178             return identity(self.getShape())
3179          else:
3180             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3181    
3182    def inverse(arg):
3183      """      """
3184      @brief      returns the inverse of the square matrix arg.
3185    
3186      @param arg      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal
3187        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3188        @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3189        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3190      """      """
3191      if isinstance(arg,escript.Data):      if isinstance(arg,numarray.NumArray):
3192         return arg.length()        return numarray.linear_algebra.inverse(arg)
3193        elif isinstance(arg,escript.Data):
3194          return escript_inverse(arg)
3195        elif isinstance(arg,float):
3196          return 1./arg
3197        elif isinstance(arg,int):
3198          return 1./float(arg)
3199        elif isinstance(arg,Symbol):
3200          return Inverse_Symbol(arg)
3201      else:      else:
3202         return sqrt((arg**2).sum())        raise TypeError,"inverse: Unknown argument type."
3203    
3204  def sign(arg):  def escript_inverse(arg): # this should be escript._inverse and use LAPACK
3205          "arg is a Data objects!!!"
3206          if not arg.getRank()==2:
3207            raise ValueError,"escript_inverse: argument must have rank 2"
3208          s=arg.getShape()      
3209          if not s[0] == s[1]:
3210            raise ValueError,"escript_inverse: argument must be a square matrix."
3211          out=escript.Data(0.,s,arg.getFunctionSpace())
3212          if s[0]==1:
3213              if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
3214                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3215              out[0,0]=1./arg[0,0]
3216          elif s[0]==2:
3217              A11=arg[0,0]
3218              A12=arg[0,1]
3219              A21=arg[1,0]
3220              A22=arg[1,1]
3221              D = A11*A22-A12*A21
3222              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3223                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3224              D=1./D
3225              out[0,0]= A22*D
3226              out[1,0]=-A21*D
3227              out[0,1]=-A12*D
3228              out[1,1]= A11*D
3229          elif s[0]==3:
3230              A11=arg[0,0]
3231              A21=arg[1,0]
3232              A31=arg[2,0]
3233              A12=arg[0,1]
3234              A22=arg[1,1]
3235              A32=arg[2,1]
3236              A13=arg[0,2]
3237              A23=arg[1,2]
3238              A33=arg[2,2]
3239              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
3240              if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
3241                  raise ZeroDivisionError,"escript_inverse: argument not invertible"
3242              D=1./D
3243              out[0,0]=(A22*A33-A23*A32)*D
3244              out[1,0]=(A31*A23-A21*A33)*D
3245              out[2,0]=(A21*A32-A31*A22)*D
3246              out[0,1]=(A13*A32-A12*A33)*D
3247              out[1,1]=(A11*A33-A31*A13)*D
3248              out[2,1]=(A12*A31-A11*A32)*D
3249              out[0,2]=(A12*A23-A13*A22)*D
3250              out[1,2]=(A13*A21-A11*A23)*D
3251              out[2,2]=(A11*A22-A12*A21)*D
3252          else:
3253             raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
3254          return out
3255    
3256    class Inverse_Symbol(DependendSymbol):
3257       """
3258       L{Symbol} representing the result of the inverse function
3259       """
3260       def __init__(self,arg):
3261          """
3262          initialization of inverse L{Symbol} with argument arg
3263          @param arg: argument of function
3264          @type arg: L{Symbol}.
3265          """
3266          if not arg.getRank()==2:
3267            raise ValueError,"Inverse_Symbol:: argument must have rank 2"
3268          s=arg.getShape()
3269          if not s[0] == s[1]:
3270            raise ValueError,"Inverse_Symbol:: argument must be a square matrix."
3271          super(Inverse_Symbol,self).__init__(args=[arg],shape=s,dim=arg.getDim())
3272