/[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 912 by gross, Wed Dec 6 03:29:49 2006 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:
# Line 238  def inf(arg): Line 239  def inf(arg):
239  #=========================================================================  #=========================================================================
240  #   some little helpers  #   some little helpers
241  #=========================================================================  #=========================================================================
242  def pokeShape(arg):  def getRank(arg):
243        """
244        identifies the rank of its argument
245    
246        @param arg: a given object
247        @type arg: L{numarray.NumArray},L{escript.Data},C{float}, C{int}, C{Symbol}
248        @return: the rank of the argument
249        @rtype: C{int}
250        @raise TypeError: if type of arg cannot be processed
251        """
252    
253        if isinstance(arg,numarray.NumArray):
254            return arg.rank
255        elif isinstance(arg,escript.Data):
256            return arg.getRank()
257        elif isinstance(arg,float):
258            return 0
259        elif isinstance(arg,int):
260            return 0
261        elif isinstance(arg,Symbol):
262            return arg.getRank()
263        else:
264          raise TypeError,"getShape: cannot identify shape"
265    def getShape(arg):
266      """      """
267      identifies the shape of its argument      identifies the shape of its argument
268    
# Line 260  def pokeShape(arg): Line 284  def pokeShape(arg):
284      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
285          return arg.getShape()          return arg.getShape()
286      else:      else:
287        raise TypeError,"pokeShape: cannot identify shape"        raise TypeError,"getShape: cannot identify shape"
288    
289  def pokeDim(arg):  def pokeDim(arg):
290      """      """
# Line 283  def commonShape(arg0,arg1): Line 307  def commonShape(arg0,arg1):
307      """      """
308      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.
309    
310      @param arg0: an object with a shape (see L{pokeShape})      @param arg0: an object with a shape (see L{getShape})
311      @param arg1: an object with a shape (see L{pokeShape})      @param arg1: an object with a shape (see L{getShape})
312      @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.
313      @rtype: C{tuple} of C{int}      @rtype: C{tuple} of C{int}
314      @raise ValueError: if no shape can be found.      @raise ValueError: if no shape can be found.
315      """      """
316      sh0=pokeShape(arg0)      sh0=getShape(arg0)
317      sh1=pokeShape(arg1)      sh1=getShape(arg1)
318      if len(sh0)<len(sh1):      if len(sh0)<len(sh1):
319         if not sh0==sh1[:len(sh0)]:         if not sh0==sh1[:len(sh0)]:
320               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 464  def matchShape(arg0,arg1):
464      @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}
465      @param arg1: a given object      @param arg1: a given object
466      @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}
467      @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.
468      @rtype: C{tuple}      @rtype: C{tuple}
469      """      """
470      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
471      sh0=pokeShape(arg0)      sh0=getShape(arg0)
472      sh1=pokeShape(arg1)      sh1=getShape(arg1)
473      if len(sh0)<len(sh):      if len(sh0)<len(sh):
474         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
475      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
# Line 573  class Symbol(object): Line 597  class Symbol(object):
597            if isinstance(a,Symbol):            if isinstance(a,Symbol):
598               out.append(a.substitute(argvals))               out.append(a.substitute(argvals))
599            else:            else:
600                s=pokeShape(s)+arg.getShape()                s=getShape(s)+arg.getShape()
601                if len(s)>0:                if len(s)>0:
602                   out.append(numarray.zeros(s),numarray.Float64)                   out.append(numarray.zeros(s),numarray.Float64)
603                else:                else:
# Line 1310  def whereNonZero(arg,tol=0.): Line 1334  def whereNonZero(arg,tol=0.):
1334     else:     else:
1335        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1336    
1337    def erf(arg):
1338       """
1339       returns erf of argument arg
1340    
1341       @param arg: argument
1342       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1343       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1344       @raises TypeError: if the type of the argument is not expected.
1345       """
1346       if isinstance(arg,escript.Data):
1347          return arg._erf()
1348       else:
1349          raise TypeError,"erf: Unknown argument type."
1350    
1351  def sin(arg):  def sin(arg):
1352     """     """
1353     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 2968  def trace(arg,axis_offset=0):
2968    
2969     @param arg: argument     @param arg: argument
2970     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2971     @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
2972                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
2973     @type axis_offset: C{int}     @type axis_offset: C{int}
2974     @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.
2975     @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 3001  def trace(arg,axis_offset=0):
3001        s=list(arg.getShape())                s=list(arg.getShape())        
3002        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3003          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)
3004        return arg._matrixtrace(axis_offset)        return arg._trace(axis_offset)
3005     elif isinstance(arg,float):     elif isinstance(arg,float):
3006        raise TypeError,"illegal argument type float."        raise TypeError,"illegal argument type float."
3007     elif isinstance(arg,int):     elif isinstance(arg,int):
# Line 2982  class Trace_Symbol(DependendSymbol): Line 3020  class Trace_Symbol(DependendSymbol):
3020        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
3021        @param arg: argument of function        @param arg: argument of function
3022        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3023        @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
3024                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3025        @type axis_offset: C{int}        @type axis_offset: C{int}
3026        """        """
3027        if arg.getRank()<2:        if arg.getRank()<2:
# Line 3049  class Trace_Symbol(DependendSymbol): Line 3087  class Trace_Symbol(DependendSymbol):
3087    
3088  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3089     """     """
3090     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.
3091    
3092     @param arg: argument     @param arg: argument
3093     @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}
3094     @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.
3095                         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.
3096     @type axis_offset: C{int}     @type axis_offset: C{int}
3097     @return: transpose of arg     @return: transpose of arg
3098     @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 3130  class Transpose_Symbol(DependendSymbol):
3130    
3131        @param arg: argument of function        @param arg: argument of function
3132        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3133        @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.
3134                         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.
3135        @type axis_offset: C{int}        @type axis_offset: C{int}
3136        """        """
3137        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3138        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3139          raise ValueError,"axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3140        s=arg.getShape()        s=arg.getShape()
3141        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())
3142    
# Line 3153  class Transpose_Symbol(DependendSymbol): Line 3191  class Transpose_Symbol(DependendSymbol):
3191           return identity(self.getShape())           return identity(self.getShape())
3192        else:        else:
3193           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3194    
3195    def swap_axes(arg,axis0=0,axis1=1):
3196       """
3197       returns the swap of arg by swaping the components axis0 and axis1
3198    
3199       @param arg: argument
3200       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3201       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3202       @type axis0: C{int}
3203       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3204       @type axis1: C{int}
3205       @return: C{arg} with swaped components
3206       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3207       """
3208       if axis0 > axis1:
3209          axis0,axis1=axis1,axis0
3210       if isinstance(arg,numarray.NumArray):
3211          return numarray.swapaxes(arg,axis0,axis1)
3212       elif isinstance(arg,escript.Data):
3213          return arg._swap_axes(axis0,axis1)
3214       elif isinstance(arg,float):
3215          raise TyepError,"float argument is not supported."
3216       elif isinstance(arg,int):
3217          raise TyepError,"int argument is not supported."
3218       elif isinstance(arg,Symbol):
3219          return SwapAxes_Symbol(arg,axis0,axis1)
3220       else:
3221          raise TypeError,"Unknown argument type."
3222    
3223    class SwapAxes_Symbol(DependendSymbol):
3224       """
3225       L{Symbol} representing the result of the swap function
3226       """
3227       def __init__(self,arg,axis0=0,axis1=1):
3228          """
3229          initialization of swap L{Symbol} with argument arg
3230    
3231          @param arg: argument
3232          @type arg: L{Symbol}.
3233          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3234          @type axis0: C{int}
3235          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3236          @type axis1: C{int}
3237          """
3238          if arg.getRank()<2:
3239             raise ValueError,"argument must have at least rank 2."
3240          if axis0<0 or axis0>arg.getRank()-1:
3241             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3242          if axis1<0 or axis1>arg.getRank()-1:
3243             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3244          if axis0 == axis1:
3245             raise ValueError,"axis indices must be different."
3246          if axis0 > axis1:
3247             axis0,axis1=axis1,axis0
3248          s=arg.getShape()
3249          s_out=[]
3250          for i in range(len(s)):
3251             if i == axis0:
3252                s_out.append(s[axis1])
3253             elif i == axis1:
3254                s_out.append(s[axis0])
3255             else:
3256                s_out.append(s[i])
3257          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3258    
3259       def getMyCode(self,argstrs,format="escript"):
3260          """
3261          returns a program code that can be used to evaluate the symbol.
3262    
3263          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3264          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3265          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3266          @type format: C{str}
3267          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3268          @rtype: C{str}
3269          @raise NotImplementedError: if the requested format is not available
3270          """
3271          if format=="escript" or format=="str"  or format=="text":
3272             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3273          else:
3274             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3275    
3276       def substitute(self,argvals):
3277          """
3278          assigns new values to symbols in the definition of the symbol.
3279          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3280    
3281          @param argvals: new values assigned to symbols
3282          @type argvals: C{dict} with keywords of type L{Symbol}.
3283          @return: result of the substitution process. Operations are executed as much as possible.
3284          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3285          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3286          """
3287          if argvals.has_key(self):
3288             arg=argvals[self]
3289             if self.isAppropriateValue(arg):
3290                return arg
3291             else:
3292                raise TypeError,"%s: new value is not appropriate."%str(self)
3293          else:
3294             arg=self.getSubstitutedArguments(argvals)
3295             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3296    
3297       def diff(self,arg):
3298          """
3299          differential of this object
3300    
3301          @param arg: the derivative is calculated with respect to arg
3302          @type arg: L{escript.Symbol}
3303          @return: derivative with respect to C{arg}
3304          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3305          """
3306          if arg==self:
3307             return identity(self.getShape())
3308          else:
3309             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3310    
3311  def symmetric(arg):  def symmetric(arg):
3312      """      """
3313      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 3683  class Add_Symbol(DependendSymbol):
3683         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3684         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3685         """         """
3686         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3687         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3688         if not sh0==sh1:         if not sh0==sh1:
3689            raise ValueError,"Add_Symbol: shape of arguments must match"            raise ValueError,"Add_Symbol: shape of arguments must match"
3690         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 3758  def mult(arg0,arg1):
3758         """         """
3759         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3760         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3761            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3762         else:         else:
3763            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3764                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3627  class Mult_Symbol(DependendSymbol): Line 3782  class Mult_Symbol(DependendSymbol):
3782         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3783         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3784         """         """
3785         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3786         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3787         if not sh0==sh1:         if not sh0==sh1:
3788            raise ValueError,"Mult_Symbol: shape of arguments must match"            raise ValueError,"Mult_Symbol: shape of arguments must match"
3789         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 3858  def quotient(arg0,arg1):
3858         """         """
3859         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3860         if testForZero(args[0]):         if testForZero(args[0]):
3861            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3862         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3863            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3864               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3732  class Quotient_Symbol(DependendSymbol): Line 3887  class Quotient_Symbol(DependendSymbol):
3887         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3888         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3889         """         """
3890         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3891         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3892         if not sh0==sh1:         if not sh0==sh1:
3893            raise ValueError,"Quotient_Symbol: shape of arguments must match"            raise ValueError,"Quotient_Symbol: shape of arguments must match"
3894         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 3964  def power(arg0,arg1):
3964         """         """
3965         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3966         if testForZero(args[0]):         if testForZero(args[0]):
3967            return numarray.zeros(pokeShape(args[0]),numarray.Float64)            return numarray.zeros(getShape(args[0]),numarray.Float64)
3968         elif testForZero(args[1]):         elif testForZero(args[1]):
3969            return numarray.ones(pokeShape(args[1]),numarray.Float64)            return numarray.ones(getShape(args[1]),numarray.Float64)
3970         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3971            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3972         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 3989  class Power_Symbol(DependendSymbol):
3989         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: if both arguments do not have the same shape.
3990         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
3991         """         """
3992         sh0=pokeShape(arg0)         sh0=getShape(arg0)
3993         sh1=pokeShape(arg1)         sh1=getShape(arg1)
3994         if not sh0==sh1:         if not sh0==sh1:
3995            raise ValueError,"Power_Symbol: shape of arguments must match"            raise ValueError,"Power_Symbol: shape of arguments must match"
3996         d0=pokeDim(arg0)         d0=pokeDim(arg0)
# Line 3934  def minimum(*args): Line 4089  def minimum(*args):
4089            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4090      return out      return out
4091    
4092  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4093      """      """
4094      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4095    
4096      @param arg: argument      @param arg: argument
4097      @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}
4098      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4099      @type minval: C{float}      @type minval: C{float} or C{None}
4100      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4101      @type maxval: C{float}      @type maxval: C{float} or C{None}
4102      @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
4103               less then maxval are unchanged.               less then maxval are unchanged.
4104      @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
4105      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4106      """      """
4107      if minval>maxval:      if not minval==None and not maxval==None:
4108         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4109      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4110        if minval == None:
4111            tmp=arg
4112        else:
4113            tmp=maximum(minval,arg)
4114        if maxval == None:
4115            return tmp
4116        else:
4117            return minimum(tmp,maxval)
4118    
4119        
4120  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3972  def inner(arg0,arg1): Line 4135  def inner(arg0,arg1):
4135      @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
4136      @raise ValueError: if the shapes of the arguments are not identical      @raise ValueError: if the shapes of the arguments are not identical
4137      """      """
4138      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4139      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4140      if not sh0==sh1:      if not sh0==sh1:
4141          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4142      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
# Line 4024  def matrix_mult(arg0,arg1): Line 4187  def matrix_mult(arg0,arg1):
4187      @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
4188      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4189      """      """
4190      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4191      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4192      if not len(sh0)==2 :      if not len(sh0)==2 :
4193          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4194      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 4236  def tensor_mult(arg0,arg1):
4236      @return: the tensor product of arg0 and arg1 at each data point      @return: the tensor product of arg0 and arg1 at each data point
4237      @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
4238      """      """
4239      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4240      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4241      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4242         return generalTensorProduct(arg0,arg1,axis_offset=1)         return generalTensorProduct(arg0,arg1,axis_offset=1)
4243      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 4311  class GeneralTensorProduct_Symbol(Depend
4311         @raise ValueError: illegal dimension         @raise ValueError: illegal dimension
4312         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4313         """         """
4314         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4315         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4316         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4317         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4318         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
# Line 4196  class GeneralTensorProduct_Symbol(Depend Line 4359  class GeneralTensorProduct_Symbol(Depend
4359           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4360           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4361    
4362  def escript_generalTensorProductNew(arg0,arg1,axis_offset):  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4363      "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!!!"
4364      # 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  
   
4365    
4366  def transposed_matrix_mult(arg0,arg1):  def transposed_matrix_mult(arg0,arg1):
4367      """      """
# Line 4275  def transposed_matrix_mult(arg0,arg1): Line 4385  def transposed_matrix_mult(arg0,arg1):
4385      @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
4386      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4387      """      """
4388      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4389      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4390      if not len(sh0)==2 :      if not len(sh0)==2 :
4391          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4392      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 4430  def transposed_tensor_mult(arg0,arg1):
4430      @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
4431      @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
4432      """      """
4433      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4434      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4435      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4436         return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)         return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4437      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 4508  class GeneralTransposedTensorProduct_Sym
4508         @raise ValueError: inconsistent dimensions of arguments.         @raise ValueError: inconsistent dimensions of arguments.
4509         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4510         """         """
4511         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4512         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4513         sh01=sh_arg0[:axis_offset]         sh01=sh_arg0[:axis_offset]
4514         sh10=sh_arg1[:axis_offset]         sh10=sh_arg1[:axis_offset]
4515         sh0=sh_arg0[axis_offset:]         sh0=sh_arg0[axis_offset:]
# Line 4448  class GeneralTransposedTensorProduct_Sym Line 4558  class GeneralTransposedTensorProduct_Sym
4558    
4559  def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct  def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4560      "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!!!"
4561      # 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  
   
4562    
4563  def matrix_transposed_mult(arg0,arg1):  def matrix_transposed_mult(arg0,arg1):
4564      """      """
# Line 4507  def matrix_transposed_mult(arg0,arg1): Line 4578  def matrix_transposed_mult(arg0,arg1):
4578      @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
4579      @raise ValueError: if the shapes of the arguments are not appropriate      @raise ValueError: if the shapes of the arguments are not appropriate
4580      """      """
4581      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4582      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4583      if not len(sh0)==2 :      if not len(sh0)==2 :
4584          raise ValueError,"first argument must have rank 2"          raise ValueError,"first argument must have rank 2"
4585      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 4614  def tensor_transposed_mult(arg0,arg1):
4614      @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
4615      @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
4616      """      """
4617      sh0=pokeShape(arg0)      sh0=getShape(arg0)
4618      sh1=pokeShape(arg1)      sh1=getShape(arg1)
4619      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):      if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4620         return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)         return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4621      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 4692  class GeneralTensorTransposedProduct_Sym
4692         @raise ValueError: inconsistent dimensions of arguments.         @raise ValueError: inconsistent dimensions of arguments.
4693         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4694         """         """
4695         sh_arg0=pokeShape(arg0)         sh_arg0=getShape(arg0)
4696         sh_arg1=pokeShape(arg1)         sh_arg1=getShape(arg1)
4697         sh0=sh_arg0[:len(sh_arg0)-axis_offset]         sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4698         sh01=sh_arg0[len(sh_arg0)-axis_offset:]         sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4699         sh10=sh_arg1[len(sh_arg1)-axis_offset:]         sh10=sh_arg1[len(sh_arg1)-axis_offset:]
# Line 4671  class GeneralTensorTransposedProduct_Sym Line 4742  class GeneralTensorTransposedProduct_Sym
4742    
4743  def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct  def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4744      "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!!!"
4745      # 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  
   
4746    
4747  #=========================================================  #=========================================================
4748  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency

Legend:
Removed from v.795  
changed lines
  Added in v.912

  ViewVC Help
Powered by ViewVC 1.1.26