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

Diff of /trunk/escript/py_src/util.py

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

revision 795 by ksteube, Mon Jul 31 01:23:58 2006 UTC revision 1042 by gross, Mon Mar 19 03:50:34 2007 UTC
# Line 26  import math Line 26  import math
26  import numarray  import numarray
27  import escript  import escript
28  import os  import os
29    from esys.escript import C_GeneralTensorProduct
30    
31  #=========================================================  #=========================================================
32  #   some helpers:  #   some helpers:
33  #=========================================================  #=========================================================
34    def getTagNames(domain):
35        """
36        returns a list of the tag names used by the domain
37    
38        
39        @param domain: a domain object
40        @type domain: C{escript.Domain}
41        @return: a list of the tag name used by the domain.
42        @rtype: C{list} of C{str}
43        """
44        return [n.strip() for n in domain.showTagNames().split(",") ]
45        
46  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
47      """      """
48      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
# Line 238  def inf(arg): Line 251  def inf(arg):
251  #=========================================================================  #=========================================================================
252  #   some little helpers  #   some little helpers
253  #=========================================================================  #=========================================================================
254  def pokeShape(arg):  def getRank(arg):
255        """
256        identifies the rank of its argument
257    
258        @param arg: a given object
259        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
260        @return: the rank of the argument
261        @rtype: C{int}
262        @raise TypeError: if type of arg cannot be processed
263        """
264    
265        if isinstance(arg,numarray.NumArray):
266            return arg.rank
267        elif isinstance(arg,escript.Data):
268            return arg.getRank()
269        elif isinstance(arg,float):
270            return 0
271        elif isinstance(arg,int):
272            return 0
273        elif isinstance(arg,Symbol):
274            return arg.getRank()
275        else:
276          raise TypeError,"getShape: cannot identify shape"
277    def getShape(arg):
278      """      """
279      identifies the shape of its argument      identifies the shape of its argument
280    
# Line 260  def pokeShape(arg): Line 296  def pokeShape(arg):
296      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
297          return arg.getShape()          return arg.getShape()
298      else:      else:
299        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
300    
301  def pokeDim(arg):  def pokeDim(arg):
302      """      """
# Line 283  def commonShape(arg0,arg1): Line 319  def commonShape(arg0,arg1):
319      """      """
320      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.      returns a shape to which arg0 can be extendent from the right and arg1 can be extended from the left.
321    
322      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
323      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
324      @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.      @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.
325      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
326      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
327      """      """
328      sh0=pokeShape(arg0)      sh0=getShape(arg0)
329      sh1=pokeShape(arg1)      sh1=getShape(arg1)
330      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
331         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
332               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"               raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
# Line 440  def matchShape(arg0,arg1): Line 476  def matchShape(arg0,arg1):
476      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}      @type arg0: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
477      @param arg1: a given object      @param arg1: a given object
478      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}      @type arg1: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, L{Symbol}
479      @return: arg0 and arg1 where copies are returned when the shape has to be changed.      @return: C{arg0} and C{arg1} where copies are returned when the shape has to be changed.
480      @rtype: C{tuple}      @rtype: C{tuple}
481      """      """
482      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
483      sh0=pokeShape(arg0)      sh0=getShape(arg0)
484      sh1=pokeShape(arg1)      sh1=getShape(arg1)
485      if len(sh0)<len(sh):      if len(sh0)<len(sh):
486         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
487      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 573  class Symbol(object): Line 609  class Symbol(object):
609            if isinstance(a,Symbol):            if isinstance(a,Symbol):
610               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
611            else:            else:
612                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
613                if len(s)>0:                if len(s)>0:
614                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
615                else:                else:
# Line 1310  def whereNonZero(arg,tol=0.): Line 1346  def whereNonZero(arg,tol=0.):
1346     else:     else:
1347        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1348    
1349    def erf(arg):
1350       """
1351       returns erf of argument arg
1352    
1353       @param arg: argument
1354       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1355       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1356       @raises TypeError: if the type of the argument is not expected.
1357       """
1358       if isinstance(arg,escript.Data):
1359          return arg._erf()
1360       else:
1361          raise TypeError,"erf: Unknown argument type."
1362    
1363  def sin(arg):  def sin(arg):
1364     """     """
1365     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 2980  def trace(arg,axis_offset=0):
2980    
2981     @param arg: argument     @param arg: argument
2982     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2983     @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     @param axis_offset: C{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
2984                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
2985     @type axis_offset: C{int}     @type axis_offset: C{int}
2986     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.     @return: trace of arg. The rank of the returned object is minus 2 of the rank of arg.
2987     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
# Line 2963  def trace(arg,axis_offset=0): Line 3013  def trace(arg,axis_offset=0):
3013        s=list(arg.getShape())                s=list(arg.getShape())        
3014        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3015          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
3016        return arg._matrixtrace(axis_offset)        return arg._trace(axis_offset)
3017     elif isinstance(arg,float):     elif isinstance(arg,float):
3018        raise TypeError,"illegal argument type float."        raise TypeError,"illegal argument type float."
3019     elif isinstance(arg,int):     elif isinstance(arg,int):
# Line 2982  class Trace_Symbol(DependendSymbol): Line 3032  class Trace_Symbol(DependendSymbol):
3032        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3033        @param arg: argument of function        @param arg: argument of function
3034        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3035        @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        @param axis_offset: C{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
3036                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3037        @type axis_offset: C{int}        @type axis_offset: C{int}
3038        """        """
3039        if arg.getRank()<2:        if arg.getRank()<2:
# Line 3049  class Trace_Symbol(DependendSymbol): Line 3099  class Trace_Symbol(DependendSymbol):
3099    
3100  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3101     """     """
3102     returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.     returns the transpose of arg by swaping the first C{axis_offset} and the last rank-axis_offset components.
3103    
3104     @param arg: argument     @param arg: argument
3105     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3106     @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.     @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3107                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3108     @type axis_offset: C{int}     @type axis_offset: C{int}
3109     @return: transpose of arg     @return: transpose of arg
3110     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.     @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
# Line 3092  class Transpose_Symbol(DependendSymbol): Line 3142  class Transpose_Symbol(DependendSymbol):
3142    
3143        @param arg: argument of function        @param arg: argument of function
3144        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3145        @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.        @param axis_offset: the first C{axis_offset} components are swapped with rest. If C{axis_offset} must be non-negative and less or equal the rank of arg.
3146                         if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.                         if C{axis_offset} is not present C{int(r/2)} where r is the rank of arg is used.
3147        @type axis_offset: C{int}        @type axis_offset: C{int}
3148        """        """
3149        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3150        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3151          raise ValueError,"axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3152        s=arg.getShape()        s=arg.getShape()
3153        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())        super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3154    
# Line 3153  class Transpose_Symbol(DependendSymbol): Line 3203  class Transpose_Symbol(DependendSymbol):
3203           return identity(self.getShape())           return identity(self.getShape())
3204        else:        else:
3205           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3206    
3207    def swap_axes(arg,axis0=0,axis1=1):
3208       """
3209       returns the swap of arg by swaping the components axis0 and axis1
3210    
3211       @param arg: argument
3212       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3213       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3214       @type axis0: C{int}
3215       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3216       @type axis1: C{int}
3217       @return: C{arg} with swaped components
3218       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3219       """
3220       if axis0 > axis1:
3221          axis0,axis1=axis1,axis0
3222       if isinstance(arg,numarray.NumArray):
3223          return numarray.swapaxes(arg,axis0,axis1)
3224       elif isinstance(arg,escript.Data):
3225          return arg._swap_axes(axis0,axis1)
3226       elif isinstance(arg,float):
3227          raise TyepError,"float argument is not supported."
3228       elif isinstance(arg,int):
3229          raise TyepError,"int argument is not supported."
3230       elif isinstance(arg,Symbol):
3231          return SwapAxes_Symbol(arg,axis0,axis1)
3232       else:
3233          raise TypeError,"Unknown argument type."
3234    
3235    class SwapAxes_Symbol(DependendSymbol):
3236       """
3237       L{Symbol} representing the result of the swap function
3238       """
3239       def __init__(self,arg,axis0=0,axis1=1):
3240          """
3241          initialization of swap L{Symbol} with argument arg
3242    
3243          @param arg: argument
3244          @type arg: L{Symbol}.
3245          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3246          @type axis0: C{int}
3247          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3248          @type axis1: C{int}
3249          """
3250          if arg.getRank()<2:
3251             raise ValueError,"argument must have at least rank 2."
3252          if axis0<0 or axis0>arg.getRank()-1:
3253             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3254          if axis1<0 or axis1>arg.getRank()-1:
3255             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3256          if axis0 == axis1:
3257             raise ValueError,"axis indices must be different."
3258          if axis0 > axis1:
3259             axis0,axis1=axis1,axis0
3260          s=arg.getShape()
3261          s_out=[]
3262          for i in range(len(s)):
3263             if i == axis0:
3264                s_out.append(s[axis1])
3265             elif i == axis1:
3266                s_out.append(s[axis0])
3267             else:
3268                s_out.append(s[i])
3269          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3270    
3271       def getMyCode(self,argstrs,format="escript"):
3272          """
3273          returns a program code that can be used to evaluate the symbol.
3274    
3275          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3276          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3277          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3278          @type format: C{str}
3279          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3280          @rtype: C{str}
3281          @raise NotImplementedError: if the requested format is not available
3282          """
3283          if format=="escript" or format=="str"  or format=="text":
3284             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3285          else:
3286             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3287    
3288       def substitute(self,argvals):
3289          """
3290          assigns new values to symbols in the definition of the symbol.
3291          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3292    
3293          @param argvals: new values assigned to symbols
3294          @type argvals: C{dict} with keywords of type L{Symbol}.
3295          @return: result of the substitution process. Operations are executed as much as possible.
3296          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3297          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3298          """
3299          if argvals.has_key(self):
3300             arg=argvals[self]
3301             if self.isAppropriateValue(arg):
3302                return arg
3303             else:
3304                raise TypeError,"%s: new value is not appropriate."%str(self)
3305          else:
3306             arg=self.getSubstitutedArguments(argvals)
3307             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3308    
3309       def diff(self,arg):
3310          """
3311          differential of this object
3312    
3313          @param arg: the derivative is calculated with respect to arg
3314          @type arg: L{escript.Symbol}
3315          @return: derivative with respect to C{arg}
3316          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3317          """
3318          if arg==self:
3319             return identity(self.getShape())
3320          else:
3321             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3322    
3323  def symmetric(arg):  def symmetric(arg):
3324      """      """
3325      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2      returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
# Line 3528  class Add_Symbol(DependendSymbol): Line 3695  class Add_Symbol(DependendSymbol):
3695         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3696         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3697         """         """
3698         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3699         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3700         if not sh0==sh1:         if not sh0==sh1:
3701            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3702         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3603  def mult(arg0,arg1): Line 3770  def mult(arg0,arg1):
3770         """         """
3771         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3772         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3773            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3774         else:         else:
3775            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3776                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3627  class Mult_Symbol(DependendSymbol): Line 3794  class Mult_Symbol(DependendSymbol):
3794         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3795         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3796         """         """
3797         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3798         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3799         if not sh0==sh1:         if not sh0==sh1:
3800            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3801         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3703  def quotient(arg0,arg1): Line 3870  def quotient(arg0,arg1):
3870         """         """
3871         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3872         if testForZero(args[0]):         if testForZero(args[0]):
3873            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3874         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3875            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3876               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3732  class Quotient_Symbol(DependendSymbol): Line 3899  class Quotient_Symbol(DependendSymbol):
3899         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3900         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3901         """         """
3902         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3903         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3904         if not sh0==sh1:         if not sh0==sh1:
3905            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3906         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])         DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0,args=[arg0,arg1])
# Line 3809  def power(arg0,arg1): Line 3976  def power(arg0,arg1):
3976         """         """
3977         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3978         if testForZero(args[0]):         if testForZero(args[0]):
3979            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3980         elif testForZero(args[1]):         elif testForZero(args[1]):
3981            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
3982         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3983            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3984         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
# Line 3834  class Power_Symbol(DependendSymbol): Line 4001  class Power_Symbol(DependendSymbol):
4001         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
4002         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4003         """         """
4004         sh0=pokeShape(arg0)         sh0=getShape(arg0)
4005         sh1=pokeShape(arg1)         sh1=getShape(arg1)
4006         if not sh0==sh1:         if not sh0==sh1:
4007            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
4008         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3934  def minimum(*args): Line 4101  def minimum(*args):
4101            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4102      return out      return out
4103    
4104  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4105      """      """
4106      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4107    
4108      @param arg: argument      @param arg: argument
4109      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float}
4110      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4111      @type minval: C{float}      @type minval: C{float} or C{None}
4112      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4113      @type maxval: C{float}      @type maxval: C{float} or C{None}
4114      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and      @return: is on object with all its value between minval and maxval. value of the argument that greater then minval and
4115               less then maxval are unchanged.               less then maxval are unchanged.
4116      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{int} or C{float} depending on the input
4117      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4118      """      """
4119      if minval>maxval:      if not minval==None and not maxval==None:
4120         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4121      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4122        if minval == None:
4123            tmp=arg
4124        else:
4125            tmp=maximum(minval,arg)
4126        if maxval == None:
4127            return tmp
4128        else:
4129            return minimum(tmp,maxval)
4130    
4131        
4132  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3972  def inner(arg0,arg1): Line 4147  def inner(arg0,arg1):
4147      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float} depending on the input
4148      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4149      """      """
4150      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4151      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4152      if not sh0==sh1:      if not sh0==sh1:
4153          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4154      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
# Line 4024  def matrix_mult(arg0,arg1): Line 4199  def matrix_mult(arg0,arg1):
4199      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4200      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4201      """      """
4202      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4203      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4204      if not len(sh0)==2 :      if not len(sh0)==2 :
4205          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4206      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
# Line 4073  def tensor_mult(arg0,arg1): Line 4248  def tensor_mult(arg0,arg1):
4248      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4249      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4250      """      """
4251      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4252      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4253      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4254         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4255      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
# Line 4148  class GeneralTensorProduct_Symbol(Depend Line 4323  class GeneralTensorProduct_Symbol(Depend
4323         @raise ValueError: illegal dimension         @raise ValueError: illegal dimension
4324         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4325         """         """
4326         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4327         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4328         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4329         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4330         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4196  class GeneralTensorProduct_Symbol(Depend Line 4371  class GeneralTensorProduct_Symbol(Depend
4371           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4372           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4373    
4374  def escript_generalTensorProductNew(arg0,arg1,axis_offset):  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4375      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4376      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  
     shape10=arg1.getShape()[:axis_offset]  
     shape1=arg1.getShape()[axis_offset:]  
     if not shape01==shape10:  
         raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)  
     # Figure out which functionspace to use (look at where operator+ is defined maybe in BinaryOp.h to get the logic for this)  
     # fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     out=GeneralTensorProduct(arg0, arg1, axis_offset)  
     return out  
   
 def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  
     "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"  
     # calculate the return shape:  
     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  
     shape10=arg1.getShape()[:axis_offset]  
     shape1=arg1.getShape()[axis_offset:]  
     if not shape01==shape10:  
         raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)  
   
     # whatr function space should be used? (this here is not good!)  
     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     # create return value:  
     out=escript.Data(0.,tuple(shape0+shape1),fs)  
     #  
     s0=[[]]  
     for k in shape0:  
           s=[]  
           for j in s0:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s0=s  
     s1=[[]]  
     for k in shape1:  
           s=[]  
           for j in s1:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s1=s  
     s01=[[]]  
     for k in shape01:  
           s=[]  
           for j in s01:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s01=s  
   
     for i0 in s0:  
        for i1 in s1:  
          s=escript.Scalar(0.,fs)  
          for i01 in s01:  
             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
   
4377    
4378  def transposed_matrix_mult(arg0,arg1):  def transposed_matrix_mult(arg0,arg1):
4379      """      """
# Line 4275  def transposed_matrix_mult(arg0,arg1): Line 4397  def transposed_matrix_mult(arg0,arg1):
4397      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4398      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4399      """      """
4400      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4401      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4402      if not len(sh0)==2 :      if not len(sh0)==2 :
4403          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4404      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
# Line 4320  def transposed_tensor_mult(arg0,arg1): Line 4442  def transposed_tensor_mult(arg0,arg1):
4442      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4443      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4444      """      """
4445      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4446      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4447      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4448         return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)         return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4449      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
# Line 4398  class GeneralTransposedTensorProduct_Sym Line 4520  class GeneralTransposedTensorProduct_Sym
4520         @raise ValueError: inconsistent dimensions of arguments.         @raise ValueError: inconsistent dimensions of arguments.
4521         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4522         """         """
4523         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4524         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4525         sh01=sh_arg0[:axis_offset]         sh01=sh_arg0[:axis_offset]
4526         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
4527         sh0=sh_arg0[axis_offset:]         sh0=sh_arg0[axis_offset:]
# Line 4448  class GeneralTransposedTensorProduct_Sym Line 4570  class GeneralTransposedTensorProduct_Sym
4570    
4571  def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct  def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4572      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4573      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
     shape01=arg0.getShape()[:axis_offset]  
     shape10=arg1.getShape()[:axis_offset]  
     shape0=arg0.getShape()[axis_offset:]  
     shape1=arg1.getShape()[axis_offset:]  
     if not shape01==shape10:  
         raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)  
   
     # whatr function space should be used? (this here is not good!)  
     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     # create return value:  
     out=escript.Data(0.,tuple(shape0+shape1),fs)  
     #  
     s0=[[]]  
     for k in shape0:  
           s=[]  
           for j in s0:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s0=s  
     s1=[[]]  
     for k in shape1:  
           s=[]  
           for j in s1:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s1=s  
     s01=[[]]  
     for k in shape01:  
           s=[]  
           for j in s01:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s01=s  
   
     for i0 in s0:  
        for i1 in s1:  
          s=escript.Scalar(0.,fs)  
          for i01 in s01:  
             s+=arg0.__getitem__(tuple(i01+i0))*arg1.__getitem__(tuple(i01+i1))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
   
4574    
4575  def matrix_transposed_mult(arg0,arg1):  def matrix_transposed_mult(arg0,arg1):
4576      """      """
# Line 4507  def matrix_transposed_mult(arg0,arg1): Line 4590  def matrix_transposed_mult(arg0,arg1):
4590      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4591      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4592      """      """
4593      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4594      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4595      if not len(sh0)==2 :      if not len(sh0)==2 :
4596          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4597      if not len(sh1)==2 and not len(sh1)==1:      if not len(sh1)==2 and not len(sh1)==1:
# Line 4543  def tensor_transposed_mult(arg0,arg1): Line 4626  def tensor_transposed_mult(arg0,arg1):
4626      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point      @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4627      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4628      """      """
4629      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4630      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4631      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4632         return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)         return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4633      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):      elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
# Line 4621  class GeneralTensorTransposedProduct_Sym Line 4704  class GeneralTensorTransposedProduct_Sym
4704         @raise ValueError: inconsistent dimensions of arguments.         @raise ValueError: inconsistent dimensions of arguments.
4705         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4706         """         """
4707         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4708         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4709         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4710         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4711         sh10=sh_arg1[len(sh_arg1)-axis_offset:]         sh10=sh_arg1[len(sh_arg1)-axis_offset:]
# Line 4671  class GeneralTensorTransposedProduct_Sym Line 4754  class GeneralTensorTransposedProduct_Sym
4754    
4755  def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct  def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4756      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"      "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4757      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
     shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  
     shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
     shape10=arg1.getShape()[arg1.getRank()-axis_offset:]  
     shape1=arg1.getShape()[:arg1.getRank()-axis_offset]  
     if not shape01==shape10:  
         raise ValueError,"dimensions of first %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)  
   
     # whatr function space should be used? (this here is not good!)  
     fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
     # create return value:  
     out=escript.Data(0.,tuple(shape0+shape1),fs)  
     #  
     s0=[[]]  
     for k in shape0:  
           s=[]  
           for j in s0:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s0=s  
     s1=[[]]  
     for k in shape1:  
           s=[]  
           for j in s1:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s1=s  
     s01=[[]]  
     for k in shape01:  
           s=[]  
           for j in s01:  
                 for i in range(k): s.append(j+[slice(i,i)])  
           s01=s  
   
     for i0 in s0:  
        for i1 in s1:  
          s=escript.Scalar(0.,fs)  
          for i01 in s01:  
             s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i1+i01))  
          out.__setitem__(tuple(i0+i1),s)  
     return out  
   
4758    
4759  #=========================================================  #=========================================================
4760  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency
Legend:
Removed from v.795  
changed lines
  Added in v.1042
  ViewVC Help
Powered by ViewVC 1.1.26