/[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 775 by ksteube, Mon Jul 10 04:00:08 2006 UTC revision 881 by gross, Thu Oct 26 02:54:47 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 440  def matchShape(arg0,arg1): Line 441  def matchShape(arg0,arg1):
441      @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}
442      @param arg1: a given object      @param arg1: a given object
443      @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}
444      @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.
445      @rtype: C{tuple}      @rtype: C{tuple}
446      """      """
447      sh=commonShape(arg0,arg1)      sh=commonShape(arg0,arg1)
# Line 1310  def whereNonZero(arg,tol=0.): Line 1311  def whereNonZero(arg,tol=0.):
1311     else:     else:
1312        raise TypeError,"whereNonZero: Unknown argument type."        raise TypeError,"whereNonZero: Unknown argument type."
1313    
1314    def erf(arg):
1315       """
1316       returns erf of argument arg
1317    
1318       @param arg: argument
1319       @type arg: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
1320       @rtype: C{float}, L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
1321       @raises TypeError: if the type of the argument is not expected.
1322       """
1323       if isinstance(arg,escript.Data):
1324          return arg._erf()
1325       else:
1326          raise TypeError,"erf: Unknown argument type."
1327    
1328  def sin(arg):  def sin(arg):
1329     """     """
1330     returns sine of argument arg     returns sine of argument arg
# Line 2930  def trace(arg,axis_offset=0): Line 2945  def trace(arg,axis_offset=0):
2945    
2946     @param arg: argument     @param arg: argument
2947     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.     @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
2948     @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
2949                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
2950     @type axis_offset: C{int}     @type axis_offset: C{int}
2951     @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.
2952     @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 2939  def trace(arg,axis_offset=0): Line 2954  def trace(arg,axis_offset=0):
2954     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
2955        sh=arg.shape        sh=arg.shape
2956        if len(sh)<2:        if len(sh)<2:
2957          raise ValueError,"trace: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
2958        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
2959          raise ValueError,"trace: axis_offset must be between 0 and %s"%len(sh)-2          raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2
2960        s1=1        s1=1
2961        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
2962        s2=1        s2=1
2963        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
2964        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2965          raise ValueError,"trace: 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)
2966        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
2967        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
2968        for i1 in range(s1):        for i1 in range(s1):
# Line 2956  def trace(arg,axis_offset=0): Line 2971  def trace(arg,axis_offset=0):
2971        out.resize(sh[:axis_offset]+sh[axis_offset+2:])        out.resize(sh[:axis_offset]+sh[axis_offset+2:])
2972        return out        return out
2973     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
2974        return escript_trace(arg,axis_offset)        if arg.getRank()<2:
2975            raise ValueError,"rank of argument must be greater than 1"
2976          if axis_offset<0 or axis_offset>arg.getRank()-2:
2977            raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
2978          s=list(arg.getShape())        
2979          if not s[axis_offset] == s[axis_offset+1]:
2980            raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2981          return arg._trace(axis_offset)
2982     elif isinstance(arg,float):     elif isinstance(arg,float):
2983        raise TypeError,"trace: illegal argument type float."        raise TypeError,"illegal argument type float."
2984     elif isinstance(arg,int):     elif isinstance(arg,int):
2985        raise TypeError,"trace: illegal argument type int."        raise TypeError,"illegal argument type int."
2986     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
2987        return Trace_Symbol(arg,axis_offset)        return Trace_Symbol(arg,axis_offset)
2988     else:     else:
2989        raise TypeError,"trace: Unknown argument type."        raise TypeError,"Unknown argument type."
2990    
 def escript_trace(arg,axis_offset):  
       "arg is a Data object"  
       if arg.getRank()<2:  
         raise ValueError,"escript_trace: rank of argument must be greater than 1"  
       if axis_offset<0 or axis_offset>arg.getRank()-2:  
         raise ValueError,"escript_trace: axis_offset must be between 0 and %s"%arg.getRank()-2  
       s=list(arg.getShape())          
       if not s[axis_offset] == s[axis_offset+1]:  
         raise ValueError,"escript_trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)  
       out=arg._matrixtrace(axis_offset)  
       return out  
2991  class Trace_Symbol(DependendSymbol):  class Trace_Symbol(DependendSymbol):
2992     """     """
2993     L{Symbol} representing the result of the trace function     L{Symbol} representing the result of the trace function
# Line 2986  class Trace_Symbol(DependendSymbol): Line 2997  class Trace_Symbol(DependendSymbol):
2997        initialization of trace L{Symbol} with argument arg        initialization of trace L{Symbol} with argument arg
2998        @param arg: argument of function        @param arg: argument of function
2999        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3000        @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
3001                    axis_offset and axis_offset+1 must be equal.                    C{axis_offset} and axis_offset+1 must be equal.
3002        @type axis_offset: C{int}        @type axis_offset: C{int}
3003        """        """
3004        if arg.getRank()<2:        if arg.getRank()<2:
3005          raise ValueError,"Trace_Symbol: rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3006        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3007          raise ValueError,"Trace_Symbol: axis_offset must be between 0 and %s"%arg.getRank()-2          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2
3008        s=list(arg.getShape())                s=list(arg.getShape())        
3009        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3010          raise ValueError,"Trace_Symbol: 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)
3011        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3012    
3013     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 3053  class Trace_Symbol(DependendSymbol): Line 3064  class Trace_Symbol(DependendSymbol):
3064    
3065  def transpose(arg,axis_offset=None):  def transpose(arg,axis_offset=None):
3066     """     """
3067     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.
3068    
3069     @param arg: argument     @param arg: argument
3070     @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}
3071     @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.
3072                         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.
3073     @type axis_offset: C{int}     @type axis_offset: C{int}
3074     @return: transpose of arg     @return: transpose of arg
3075     @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 3067  def transpose(arg,axis_offset=None): Line 3078  def transpose(arg,axis_offset=None):
3078        if axis_offset==None: axis_offset=int(arg.rank/2)        if axis_offset==None: axis_offset=int(arg.rank/2)
3079        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))        return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3080     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
3081        if axis_offset==None: axis_offset=int(arg.getRank()/2)        r=arg.getRank()
3082        return escript_transpose(arg,axis_offset)        if axis_offset==None: axis_offset=int(r/2)
3083          if axis_offset<0 or axis_offset>r:
3084            raise ValueError,"axis_offset must be between 0 and %s"%r
3085          return arg._transpose(axis_offset)
3086     elif isinstance(arg,float):     elif isinstance(arg,float):
3087        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3088          raise ValueError,"transpose: axis_offset must be 0 for float argument"          raise ValueError,"axis_offset must be 0 for float argument"
3089        return arg        return arg
3090     elif isinstance(arg,int):     elif isinstance(arg,int):
3091        if not ( axis_offset==0 or axis_offset==None):        if not ( axis_offset==0 or axis_offset==None):
3092          raise ValueError,"transpose: axis_offset must be 0 for int argument"          raise ValueError,"axis_offset must be 0 for int argument"
3093        return float(arg)        return float(arg)
3094     elif isinstance(arg,Symbol):     elif isinstance(arg,Symbol):
3095        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3096        return Transpose_Symbol(arg,axis_offset)        return Transpose_Symbol(arg,axis_offset)
3097     else:     else:
3098        raise TypeError,"transpose: Unknown argument type."        raise TypeError,"Unknown argument type."
3099    
 def escript_transpose(arg,axis_offset):  
       "arg is a Data object"  
       r=arg.getRank()  
       if axis_offset<0 or axis_offset>r:  
         raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r  
       out=arg._transpose(axis_offset)  
       return out  
3100  class Transpose_Symbol(DependendSymbol):  class Transpose_Symbol(DependendSymbol):
3101     """     """
3102     L{Symbol} representing the result of the transpose function     L{Symbol} representing the result of the transpose function
# Line 3100  class Transpose_Symbol(DependendSymbol): Line 3107  class Transpose_Symbol(DependendSymbol):
3107    
3108        @param arg: argument of function        @param arg: argument of function
3109        @type arg: L{Symbol}.        @type arg: L{Symbol}.
3110        @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.
3111                         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.
3112        @type axis_offset: C{int}        @type axis_offset: C{int}
3113        """        """
3114        if axis_offset==None: axis_offset=int(arg.getRank()/2)        if axis_offset==None: axis_offset=int(arg.getRank()/2)
3115        if axis_offset<0 or axis_offset>arg.getRank():        if axis_offset<0 or axis_offset>arg.getRank():
3116          raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()
3117        s=arg.getShape()        s=arg.getShape()
3118        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())
3119    
# Line 3161  class Transpose_Symbol(DependendSymbol): Line 3168  class Transpose_Symbol(DependendSymbol):
3168           return identity(self.getShape())           return identity(self.getShape())
3169        else:        else:
3170           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3171    
3172    def swap_axes(arg,axis0=0,axis1=1):
3173       """
3174       returns the swap of arg by swaping the components axis0 and axis1
3175    
3176       @param arg: argument
3177       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}.
3178       @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3179       @type axis0: C{int}
3180       @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3181       @type axis1: C{int}
3182       @return: C{arg} with swaped components
3183       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray} depending on the type of arg.
3184       """
3185       if axis0 > axis1:
3186          axis0,axis1=axis1,axis0
3187       if isinstance(arg,numarray.NumArray):
3188          return numarray.swapaxes(arg,axis0,axis1)
3189       elif isinstance(arg,escript.Data):
3190          return arg._swap_axes(axis0,axis1)
3191       elif isinstance(arg,float):
3192          raise TyepError,"float argument is not supported."
3193       elif isinstance(arg,int):
3194          raise TyepError,"int argument is not supported."
3195       elif isinstance(arg,Symbol):
3196          return SwapAxes_Symbol(arg,axis0,axis1)
3197       else:
3198          raise TypeError,"Unknown argument type."
3199    
3200    class SwapAxes_Symbol(DependendSymbol):
3201       """
3202       L{Symbol} representing the result of the swap function
3203       """
3204       def __init__(self,arg,axis0=0,axis1=1):
3205          """
3206          initialization of swap L{Symbol} with argument arg
3207    
3208          @param arg: argument
3209          @type arg: L{Symbol}.
3210          @param axis0: axis. C{axis0} must be non-negative and less than the rank of arg.
3211          @type axis0: C{int}
3212          @param axis1: axis. C{axis1} must be non-negative and less than the rank of arg.
3213          @type axis1: C{int}
3214          """
3215          if arg.getRank()<2:
3216             raise ValueError,"argument must have at least rank 2."
3217          if axis0<0 or axis0>arg.getRank()-1:
3218             raise ValueError,"axis0 must be between 0 and %s"%arg.getRank()-1
3219          if axis1<0 or axis1>arg.getRank()-1:
3220             raise ValueError,"axis1 must be between 0 and %s"%arg.getRank()-1
3221          if axis0 == axis1:
3222             raise ValueError,"axis indices must be different."
3223          if axis0 > axis1:
3224             axis0,axis1=axis1,axis0
3225          s=arg.getShape()
3226          s_out=[]
3227          for i in range(len(s)):
3228             if i == axis0:
3229                s_out.append(s[axis1])
3230             elif i == axis1:
3231                s_out.append(s[axis0])
3232             else:
3233                s_out.append(s[i])
3234          super(SwapAxes_Symbol,self).__init__(args=[arg,axis0,axis1],shape=tuple(s_out),dim=arg.getDim())
3235    
3236       def getMyCode(self,argstrs,format="escript"):
3237          """
3238          returns a program code that can be used to evaluate the symbol.
3239    
3240          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3241          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3242          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3243          @type format: C{str}
3244          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3245          @rtype: C{str}
3246          @raise NotImplementedError: if the requested format is not available
3247          """
3248          if format=="escript" or format=="str"  or format=="text":
3249             return "swap(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3250          else:
3251             raise NotImplementedError,"SwapAxes_Symbol does not provide program code for format %s."%format
3252    
3253       def substitute(self,argvals):
3254          """
3255          assigns new values to symbols in the definition of the symbol.
3256          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3257    
3258          @param argvals: new values assigned to symbols
3259          @type argvals: C{dict} with keywords of type L{Symbol}.
3260          @return: result of the substitution process. Operations are executed as much as possible.
3261          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3262          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3263          """
3264          if argvals.has_key(self):
3265             arg=argvals[self]
3266             if self.isAppropriateValue(arg):
3267                return arg
3268             else:
3269                raise TypeError,"%s: new value is not appropriate."%str(self)
3270          else:
3271             arg=self.getSubstitutedArguments(argvals)
3272             return swap_axes(arg[0],axis0=arg[1],axis1=arg[2])
3273    
3274       def diff(self,arg):
3275          """
3276          differential of this object
3277    
3278          @param arg: the derivative is calculated with respect to arg
3279          @type arg: L{escript.Symbol}
3280          @return: derivative with respect to C{arg}
3281          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3282          """
3283          if arg==self:
3284             return identity(self.getShape())
3285          else:
3286             return swap_axes(self.getDifferentiatedArguments(arg)[0],axis0=self.getArgument()[1],axis1=self.getArgument()[2])
3287    
3288  def symmetric(arg):  def symmetric(arg):
3289      """      """
3290      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 3173  def symmetric(arg): Line 3297  def symmetric(arg):
3297      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3298        if arg.rank==2:        if arg.rank==2:
3299          if not (arg.shape[0]==arg.shape[1]):          if not (arg.shape[0]==arg.shape[1]):
3300             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3301        elif arg.rank==4:        elif arg.rank==4:
3302          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):          if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3303             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3304        else:        else:
3305          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3306        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3307      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3308        return escript_symmetric(arg)        if arg.getRank()==2:
3309            if not (arg.getShape()[0]==arg.getShape()[1]):
3310               raise ValueError,"argument must be square."
3311            return arg._symmetric()
3312          elif arg.getRank()==4:
3313            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3314               raise ValueError,"argument must be square."
3315            return arg._symmetric()
3316          else:
3317            raise ValueError,"rank 2 or 4 is required."
3318      elif isinstance(arg,float):      elif isinstance(arg,float):
3319        return arg        return arg
3320      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3189  def symmetric(arg): Line 3322  def symmetric(arg):
3322      elif isinstance(arg,Symbol):      elif isinstance(arg,Symbol):
3323        if arg.getRank()==2:        if arg.getRank()==2:
3324          if not (arg.getShape()[0]==arg.getShape()[1]):          if not (arg.getShape()[0]==arg.getShape()[1]):
3325             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3326        elif arg.getRank()==4:        elif arg.getRank()==4:
3327          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):          if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3328             raise ValueError,"symmetric: argument must be square."             raise ValueError,"argument must be square."
3329        else:        else:
3330          raise ValueError,"symmetric: rank 2 or 4 is required."          raise ValueError,"rank 2 or 4 is required."
3331        return (arg+transpose(arg))/2        return (arg+transpose(arg))/2
3332      else:      else:
3333        raise TypeError,"symmetric: Unknown argument type."        raise TypeError,"symmetric: Unknown argument type."
3334    
 def escript_symmetric(arg):  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=arg._symmetric()  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_symmetric: argument must be square."  
         out=arg._symmetric()  
       else:  
         raise ValueError,"escript_symmetric: rank 2 or 4 is required."  
       return out  
   
3335  def nonsymmetric(arg):  def nonsymmetric(arg):
3336      """      """
3337      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2      returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
# Line 3232  def nonsymmetric(arg): Line 3352  def nonsymmetric(arg):
3352          raise ValueError,"nonsymmetric: rank 2 or 4 is required."          raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3353        return (arg-transpose(arg))/2        return (arg-transpose(arg))/2
3354      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
3355        return escript_nonsymmetric(arg)        if arg.getRank()==2:
3356            if not (arg.getShape()[0]==arg.getShape()[1]):
3357               raise ValueError,"argument must be square."
3358            return arg._nonsymmetric()
3359          elif arg.getRank()==4:
3360            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3361               raise ValueError,"argument must be square."
3362            return arg._nonsymmetric()
3363          else:
3364            raise ValueError,"rank 2 or 4 is required."
3365      elif isinstance(arg,float):      elif isinstance(arg,float):
3366        return arg        return arg
3367      elif isinstance(arg,int):      elif isinstance(arg,int):
# Line 3250  def nonsymmetric(arg): Line 3379  def nonsymmetric(arg):
3379      else:      else:
3380        raise TypeError,"nonsymmetric: Unknown argument type."        raise TypeError,"nonsymmetric: Unknown argument type."
3381    
 def escript_nonsymmetric(arg): # this should be implemented in c++  
       if arg.getRank()==2:  
         if not (arg.getShape()[0]==arg.getShape()[1]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=arg._nonsymmetric()  
       elif arg.getRank()==4:  
         if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):  
            raise ValueError,"escript_nonsymmetric: argument must be square."  
         out=arg._nonsymmetric()  
       else:  
         raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."  
       return out  
   
   
3382  def inverse(arg):  def inverse(arg):
3383      """      """
3384      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3385    
3386      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.      @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3387      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3388      @return: inverse arg_inv of the argument. It will be matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])      @return: inverse arg_inv of the argument. It will be matrix_mult(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3389      @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
3390      @note: for L{escript.Data} objects the dimension is restricted to 3.      @note: for L{escript.Data} objects the dimension is restricted to 3.
3391      """      """
# Line 3407  class Inverse_Symbol(DependendSymbol): Line 3522  class Inverse_Symbol(DependendSymbol):
3522        if arg==self:        if arg==self:
3523           return identity(self.getShape())           return identity(self.getShape())
3524        else:        else:
3525           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrix_mult(matrix_mult(self,self.getDifferentiatedArguments(arg)[0]),self)
3526    
3527  def eigenvalues(arg):  def eigenvalues(arg):
3528      """      """
# Line 3951  def minimum(*args): Line 4066  def minimum(*args):
4066            out=add(out,mult(whereNegative(diff),diff))            out=add(out,mult(whereNegative(diff),diff))
4067      return out      return out
4068    
4069  def clip(arg,minval=0.,maxval=1.):  def clip(arg,minval=None,maxval=None):
4070      """      """
4071      cuts the values of arg between minval and maxval      cuts the values of arg between minval and maxval
4072    
4073      @param arg: argument      @param arg: argument
4074      @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}
4075      @param minval: lower range      @param minval: lower range. If None no lower range is applied
4076      @type minval: C{float}      @type minval: C{float} or C{None}
4077      @param maxval: upper range      @param maxval: upper range. If None no upper range is applied
4078      @type maxval: C{float}      @type maxval: C{float} or C{None}
4079      @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
4080               less then maxval are unchanged.               less then maxval are unchanged.
4081      @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
4082      @raise ValueError: if minval>maxval      @raise ValueError: if minval>maxval
4083      """      """
4084      if minval>maxval:      if not minval==None and not maxval==None:
4085         raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)         if minval>maxval:
4086      return minimum(maximum(minval,arg),maxval)            raise ValueError,"minval = %s must be less then maxval %s"%(minval,maxval)
4087        if minval == None:
4088            tmp=arg
4089        else:
4090            tmp=maximum(minval,arg)
4091        if maxval == None:
4092            return tmp
4093        else:
4094            return minimum(tmp,maxval)
4095    
4096        
4097  def inner(arg0,arg1):  def inner(arg0,arg1):
# Line 3995  def inner(arg0,arg1): Line 4118  def inner(arg0,arg1):
4118          raise ValueError,"inner: shape of arguments does not match"          raise ValueError,"inner: shape of arguments does not match"
4119      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))      return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
4120    
4121    def outer(arg0,arg1):
4122        """
4123        the outer product of the two argument:
4124    
4125        out[t,s]=arg0[t]*arg1[s]
4126    
4127        where
4128    
4129            - s runs through arg0.Shape
4130            - t runs through arg1.Shape
4131    
4132        @param arg0: first argument
4133        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4134        @param arg1: second argument
4135        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4136        @return: the outer product of arg0 and arg1 at each data point
4137        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4138        """
4139        return generalTensorProduct(arg0,arg1,axis_offset=0)
4140    
4141  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
4142      """      """
4143        see L{matrix_mult}
4144        """
4145        return matrix_mult(arg0,arg1)
4146    
4147    def matrix_mult(arg0,arg1):
4148        """
4149      matrix-matrix or matrix-vector product of the two argument:      matrix-matrix or matrix-vector product of the two argument:
4150    
4151      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]      out[s0]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0]
# Line 4005  def matrixmult(arg0,arg1): Line 4154  def matrixmult(arg0,arg1):
4154    
4155      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]      out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[r0,s1]
4156    
4157      The second dimension of arg0 and the length of arg1 must match.      The second dimension of arg0 and the first dimension of arg1 must match.
4158    
4159      @param arg0: first argument of rank 2      @param arg0: first argument of rank 2
4160      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4023  def matrixmult(arg0,arg1): Line 4172  def matrixmult(arg0,arg1):
4172          raise ValueError,"second argument must have rank 1 or 2"          raise ValueError,"second argument must have rank 1 or 2"
4173      return generalTensorProduct(arg0,arg1,axis_offset=1)      return generalTensorProduct(arg0,arg1,axis_offset=1)
4174    
4175  def outer(arg0,arg1):  def tensormult(arg0,arg1):
4176      """      """
4177      the outer product of the two argument:      use L{tensor_mult}
   
     out[t,s]=arg0[t]*arg1[s]  
   
     where  
   
         - s runs through arg0.Shape  
         - t runs through arg1.Shape  
   
     @param arg0: first argument  
     @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @param arg1: second argument  
     @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}  
     @return: the outer product of arg0 and arg1 at each data point  
     @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input  
4178      """      """
4179      return generalTensorProduct(arg0,arg1,axis_offset=0)      return tensor_mult(arg0,arg1)
   
4180    
4181  def tensormult(arg0,arg1):  def tensor_mult(arg0,arg1):
4182      """      """
4183      the tensor product of the two argument:      the tensor product of the two argument:
4184            
# Line 4069  def tensormult(arg0,arg1): Line 4203  def tensormult(arg0,arg1):
4203    
4204      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]      out[s0,s1]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]
4205    
4206      In the first case the the second dimension of arg0 and the length of arg1 must match and        In the first case the the second dimension of arg0 and the last dimension of arg1 must match and  
4207      in the second case the two last dimensions of arg0 must match the shape of arg1.      in the second case the two last dimensions of arg0 must match the two first dimensions of arg1.
4208    
4209      @param arg0: first argument of rank 2 or 4      @param arg0: first argument of rank 2 or 4
4210      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
# Line 4086  def tensormult(arg0,arg1): Line 4220  def tensormult(arg0,arg1):
4220      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):
4221         return generalTensorProduct(arg0,arg1,axis_offset=2)         return generalTensorProduct(arg0,arg1,axis_offset=2)
4222      else:      else:
4223          raise ValueError,"tensormult: first argument must have rank 2 or 4"          raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
4224    
4225  def generalTensorProduct(arg0,arg1,axis_offset=0):  def generalTensorProduct(arg0,arg1,axis_offset=0):
4226      """      """
# Line 4100  def generalTensorProduct(arg0,arg1,axis_ Line 4234  def generalTensorProduct(arg0,arg1,axis_
4234          - r runs trough arg0.Shape[:axis_offset]          - r runs trough arg0.Shape[:axis_offset]
4235          - t runs through arg1.Shape[axis_offset:]          - t runs through arg1.Shape[axis_offset:]
4236    
     In the first case the the second dimension of arg0 and the length of arg1 must match and    
     in the second case the two last dimensions of arg0 must match the shape of arg1.  
   
4237      @param arg0: first argument      @param arg0: first argument
4238      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4239      @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0      @param arg1: second argument
4240      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4241      @return: the general tensor product of arg0 and arg1 at each data point.      @return: the general tensor product of arg0 and arg1 at each data point.
4242      @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
# Line 4118  def generalTensorProduct(arg0,arg1,axis_ Line 4249  def generalTensorProduct(arg0,arg1,axis_
4249             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)             return GeneralTensorProduct_Symbol(arg0,arg1,axis_offset)
4250         else:         else:
4251             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:             if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[:axis_offset]:
4252                 raise ValueError,"generalTensorProduct: dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)                 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)
4253             arg0_c=arg0.copy()             arg0_c=arg0.copy()
4254             arg1_c=arg1.copy()             arg1_c=arg1.copy()
4255             sh0,sh1=arg0.shape,arg1.shape             sh0,sh1=arg0.shape,arg1.shape
# Line 4144  def generalTensorProduct(arg0,arg1,axis_ Line 4275  def generalTensorProduct(arg0,arg1,axis_
4275                                    
4276  class GeneralTensorProduct_Symbol(DependendSymbol):  class GeneralTensorProduct_Symbol(DependendSymbol):
4277     """     """
4278     Symbol representing the quotient of two arguments.     Symbol representing the general tensor product of two arguments
4279     """     """
4280     def __init__(self,arg0,arg1,axis_offset=0):     def __init__(self,arg0,arg1,axis_offset=0):
4281         """         """
4282         initialization of L{Symbol} representing the quotient of two arguments         initialization of L{Symbol} representing the general tensor product of two arguments.
4283    
4284         @param arg0: numerator         @param arg0: first argument
4285         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4286         @param arg1: denominator         @param arg1: second argument
4287         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.         @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4288         @raise ValueError: if both arguments do not have the same shape.         @raise ValueError: illegal dimension
4289         @note: if both arguments have a spatial dimension, they must equal.         @note: if both arguments have a spatial dimension, they must equal.
4290         """         """
4291         sh_arg0=pokeShape(arg0)         sh_arg0=pokeShape(arg0)
# Line 4205  class GeneralTensorProduct_Symbol(Depend Line 4336  class GeneralTensorProduct_Symbol(Depend
4336           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
4337           return generalTensorProduct(args[0],args[1],args[2])           return generalTensorProduct(args[0],args[1],args[2])
4338    
4339  def escript_generalTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorProduct  def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
4340      "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!!!"
4341      # calculate the return shape:      return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
4342      shape0=arg0.getShape()[:arg0.getRank()-axis_offset]  
4343      shape01=arg0.getShape()[arg0.getRank()-axis_offset:]  def transposed_matrix_mult(arg0,arg1):
4344      shape10=arg1.getShape()[:axis_offset]      """
4345      shape1=arg1.getShape()[axis_offset:]      transposed(matrix)-matrix or transposed(matrix)-vector product of the two argument:
4346      if not shape01==shape10:  
4347          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)      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]
4348    
4349      # whatr function space should be used? (this here is not good!)      or
4350      fs=(escript.Scalar(0.,arg0.getFunctionSpace())+escript.Scalar(0.,arg1.getFunctionSpace())).getFunctionSpace()  
4351      # create return value:      out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4352      out=escript.Data(0.,tuple(shape0+shape1),fs)  
4353      #      The function call transposed_matrix_mult(arg0,arg1) is equivalent to matrix_mult(transpose(arg0),arg1).
4354      s0=[[]]  
4355      for k in shape0:      The first dimension of arg0 and arg1 must match.
4356            s=[]  
4357            for j in s0:      @param arg0: first argument of rank 2
4358                  for i in range(k): s.append(j+[slice(i,i)])      @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4359            s0=s      @param arg1: second argument of at least rank 1
4360      s1=[[]]      @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4361      for k in shape1:      @return: the product of the transposed of arg0 and arg1 at each data point
4362            s=[]      @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4363            for j in s1:      @raise ValueError: if the shapes of the arguments are not appropriate
4364                  for i in range(k): s.append(j+[slice(i,i)])      """
4365            s1=s      sh0=pokeShape(arg0)
4366      s01=[[]]      sh1=pokeShape(arg1)
4367      for k in shape01:      if not len(sh0)==2 :
4368            s=[]          raise ValueError,"first argument must have rank 2"
4369            for j in s01:      if not len(sh1)==2 and not len(sh1)==1:
4370                  for i in range(k): s.append(j+[slice(i,i)])          raise ValueError,"second argument must have rank 1 or 2"
4371            s01=s      return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4372    
4373      for i0 in s0:  def transposed_tensor_mult(arg0,arg1):
4374         for i1 in s1:      """
4375           s=escript.Scalar(0.,fs)      the tensor product of the transposed of the first and the second argument
4376           for i01 in s01:      
4377              s+=arg0.__getitem__(tuple(i0+i01))*arg1.__getitem__(tuple(i01+i1))      for arg0 of rank 2 this is
4378           out.__setitem__(tuple(i0+i1),s)  
4379      return out      out[s0]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0]  
4380    
4381        or
4382    
4383        out[s0,s1]=S{Sigma}_{r0} arg0[r0,s0]*arg1[r0,s1]
4384    
4385      
4386        and for arg0 of rank 4 this is
4387    
4388        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]
4389                  
4390        or
4391    
4392        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]
4393    
4394        or
4395    
4396        out[s0,s1]=S{Sigma}_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]
4397    
4398        In the first case the the first dimension of arg0 and the first dimension of arg1 must match and  
4399        in the second case the two first dimensions of arg0 must match the two first dimension of arg1.
4400    
4401        The function call transposed_tensor_mult(arg0,arg1) is equivalent to tensor_mult(transpose(arg0),arg1).
4402    
4403        @param arg0: first argument of rank 2 or 4
4404        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4405        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4406        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4407        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4408        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4409        """
4410        sh0=pokeShape(arg0)
4411        sh1=pokeShape(arg1)
4412        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4413           return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
4414        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4415           return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
4416        else:
4417            raise ValueError,"first argument must have rank 2 or 4"
4418    
4419    def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
4420        """
4421        generalized tensor product of transposed of arg0 and arg1:
4422    
4423        out[s,t]=S{Sigma}_r arg0[r,s]*arg1[r,t]
4424    
4425        where
4426    
4427            - s runs through arg0.Shape[axis_offset:]
4428            - r runs trough arg0.Shape[:axis_offset]
4429            - t runs through arg1.Shape[axis_offset:]
4430    
4431        The function call generalTransposedTensorProduct(arg0,arg1,axis_offset) is equivalent
4432        to generalTensorProduct(transpose(arg0,arg0.rank-axis_offset),arg1,axis_offset).
4433    
4434        @param arg0: first argument
4435        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4436        @param arg1: second argument
4437        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4438        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4439        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4440        """
4441        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4442        arg0,arg1=matchType(arg0,arg1)
4443        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4444        if isinstance(arg0,numarray.NumArray):
4445           if isinstance(arg1,Symbol):
4446               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4447           else:
4448               if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
4449                   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)
4450               arg0_c=arg0.copy()
4451               arg1_c=arg1.copy()
4452               sh0,sh1=arg0.shape,arg1.shape
4453               d0,d1,d01=1,1,1
4454               for i in sh0[axis_offset:]: d0*=i
4455               for i in sh1[axis_offset:]: d1*=i
4456               for i in sh0[:axis_offset]: d01*=i
4457               arg0_c.resize((d01,d0))
4458               arg1_c.resize((d01,d1))
4459               out=numarray.zeros((d0,d1),numarray.Float64)
4460               for i0 in range(d0):
4461                        for i1 in range(d1):
4462                             out[i0,i1]=numarray.sum(arg0_c[:,i0]*arg1_c[:,i1])
4463               out.resize(sh0[axis_offset:]+sh1[axis_offset:])
4464               return out
4465        elif isinstance(arg0,escript.Data):
4466           if isinstance(arg1,Symbol):
4467               return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4468           else:
4469               return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4470        else:      
4471           return GeneralTransposedTensorProduct_Symbol(arg0,arg1,axis_offset)
4472                    
4473    class GeneralTransposedTensorProduct_Symbol(DependendSymbol):
4474       """
4475       Symbol representing the general tensor product of the transposed of arg0 and arg1
4476       """
4477       def __init__(self,arg0,arg1,axis_offset=0):
4478           """
4479           initialization of L{Symbol} representing tensor product of the transposed of arg0 and arg1
4480    
4481           @param arg0: first argument
4482           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4483           @param arg1: second argument
4484           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4485           @raise ValueError: inconsistent dimensions of arguments.
4486           @note: if both arguments have a spatial dimension, they must equal.
4487           """
4488           sh_arg0=pokeShape(arg0)
4489           sh_arg1=pokeShape(arg1)
4490           sh01=sh_arg0[:axis_offset]
4491           sh10=sh_arg1[:axis_offset]
4492           sh0=sh_arg0[axis_offset:]
4493           sh1=sh_arg1[axis_offset:]
4494           if not sh01==sh10:
4495               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)
4496           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4497    
4498       def getMyCode(self,argstrs,format="escript"):
4499          """
4500          returns a program code that can be used to evaluate the symbol.
4501    
4502          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4503          @type argstrs: C{list} of length 2 of C{str}.
4504          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4505          @type format: C{str}
4506          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4507          @rtype: C{str}
4508          @raise NotImplementedError: if the requested format is not available
4509          """
4510          if format=="escript" or format=="str" or format=="text":
4511             return "generalTransposedTensorProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4512          else:
4513             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4514    
4515       def substitute(self,argvals):
4516          """
4517          assigns new values to symbols in the definition of the symbol.
4518          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4519    
4520          @param argvals: new values assigned to symbols
4521          @type argvals: C{dict} with keywords of type L{Symbol}.
4522          @return: result of the substitution process. Operations are executed as much as possible.
4523          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4524          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4525          """
4526          if argvals.has_key(self):
4527             arg=argvals[self]
4528             if self.isAppropriateValue(arg):
4529                return arg
4530             else:
4531                raise TypeError,"%s: new value is not appropriate."%str(self)
4532          else:
4533             args=self.getSubstitutedArguments(argvals)
4534             return generalTransposedTensorProduct(args[0],args[1],args[2])
4535    
4536    def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset): # this should be escript._generalTransposedTensorProduct
4537        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4538        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
4539    
4540    def matrix_transposed_mult(arg0,arg1):
4541        """
4542        matrix-transposed(matrix) product of the two argument:
4543    
4544        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4545    
4546        The function call matrix_transposed_mult(arg0,arg1) is equivalent to matrix_mult(arg0,transpose(arg1)).
4547    
4548        The last dimensions of arg0 and arg1 must match.
4549    
4550        @param arg0: first argument of rank 2
4551        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4552        @param arg1: second argument of rank 2
4553        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4554        @return: the product of arg0 and the transposed of arg1 at each data point
4555        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4556        @raise ValueError: if the shapes of the arguments are not appropriate
4557        """
4558        sh0=pokeShape(arg0)
4559        sh1=pokeShape(arg1)
4560        if not len(sh0)==2 :
4561            raise ValueError,"first argument must have rank 2"
4562        if not len(sh1)==2 and not len(sh1)==1:
4563            raise ValueError,"second argument must have rank 1 or 2"
4564        return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4565    
4566    def tensor_transposed_mult(arg0,arg1):
4567        """
4568        the tensor product of the first and the transpose of the second argument
4569        
4570        for arg0 of rank 2 this is
4571    
4572        out[s0,s1]=S{Sigma}_{r0} arg0[s0,r0]*arg1[s1,r0]
4573    
4574        and for arg0 of rank 4 this is
4575    
4576        out[s0,s1,s2,s3]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]
4577                  
4578        or
4579    
4580        out[s0,s1,s2]=S{Sigma}_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]
4581    
4582        In the first case the the second dimension of arg0 and arg1 must match and  
4583        in the second case the two last dimensions of arg0 must match the two last dimension of arg1.
4584    
4585        The function call tensor_transpose_mult(arg0,arg1) is equivalent to tensor_mult(arg0,transpose(arg1)).
4586    
4587        @param arg0: first argument of rank 2 or 4
4588        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4589        @param arg1: second argument of shape greater of 1 or 2 depending on rank of arg0
4590        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
4591        @return: the tensor product of tarnsposed of arg0 and arg1 at each data point
4592        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4593        """
4594        sh0=pokeShape(arg0)
4595        sh1=pokeShape(arg1)
4596        if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
4597           return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
4598        elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
4599           return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
4600        else:
4601            raise ValueError,"first argument must have rank 2 or 4"
4602    
4603    def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
4604        """
4605        generalized tensor product of transposed of arg0 and arg1:
4606    
4607        out[s,t]=S{Sigma}_r arg0[s,r]*arg1[t,r]
4608    
4609        where
4610    
4611            - s runs through arg0.Shape[:arg0.Rank-axis_offset]
4612            - r runs trough arg0.Shape[arg1.Rank-axis_offset:]
4613            - t runs through arg1.Shape[arg1.Rank-axis_offset:]
4614    
4615        The function call generalTensorTransposedProduct(arg0,arg1,axis_offset) is equivalent
4616        to generalTensorProduct(arg0,transpose(arg1,arg1.Rank-axis_offset),axis_offset).
4617    
4618        @param arg0: first argument
4619        @type arg0: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4620        @param arg1: second argument
4621        @type arg1: L{numarray.NumArray}, L{escript.Data}, L{Symbol}, C{float}, C{int}
4622        @return: the general tensor product of transposed(arg0) and arg1 at each data point.
4623        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
4624        """
4625        if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
4626        arg0,arg1=matchType(arg0,arg1)
4627        # at this stage arg0 and arg0 are both numarray.NumArray or escript.Data or Symbols
4628        if isinstance(arg0,numarray.NumArray):
4629           if isinstance(arg1,Symbol):
4630               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4631           else:
4632               if not arg0.shape[arg0.rank-axis_offset:]==arg1.shape[arg1.rank-axis_offset:]:
4633                   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)
4634               arg0_c=arg0.copy()
4635               arg1_c=arg1.copy()
4636               sh0,sh1=arg0.shape,arg1.shape
4637               d0,d1,d01=1,1,1
4638               for i in sh0[:arg0.rank-axis_offset]: d0*=i
4639               for i in sh1[:arg1.rank-axis_offset]: d1*=i
4640               for i in sh1[arg1.rank-axis_offset:]: d01*=i
4641               arg0_c.resize((d0,d01))
4642               arg1_c.resize((d1,d01))
4643               out=numarray.zeros((d0,d1),numarray.Float64)
4644               for i0 in range(d0):
4645                        for i1 in range(d1):
4646                             out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[i1,:])
4647               out.resize(sh0[:arg0.rank-axis_offset]+sh1[:arg1.rank-axis_offset])
4648               return out
4649        elif isinstance(arg0,escript.Data):
4650           if isinstance(arg1,Symbol):
4651               return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4652           else:
4653               return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
4654        else:      
4655           return GeneralTensorTransposedProduct_Symbol(arg0,arg1,axis_offset)
4656                    
4657    class GeneralTensorTransposedProduct_Symbol(DependendSymbol):
4658       """
4659       Symbol representing the general tensor product of arg0 and the transpose of arg1
4660       """
4661       def __init__(self,arg0,arg1,axis_offset=0):
4662           """
4663           initialization of L{Symbol} representing the general tensor product of arg0 and the transpose of arg1
4664    
4665           @param arg0: first argument
4666           @type arg0: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4667           @param arg1: second argument
4668           @type arg1: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray}.
4669           @raise ValueError: inconsistent dimensions of arguments.
4670           @note: if both arguments have a spatial dimension, they must equal.
4671           """
4672           sh_arg0=pokeShape(arg0)
4673           sh_arg1=pokeShape(arg1)
4674           sh0=sh_arg0[:len(sh_arg0)-axis_offset]
4675           sh01=sh_arg0[len(sh_arg0)-axis_offset:]
4676           sh10=sh_arg1[len(sh_arg1)-axis_offset:]
4677           sh1=sh_arg1[:len(sh_arg1)-axis_offset]
4678           if not sh01==sh10:
4679               raise ValueError,"dimensions of last %s components in left argument don't match the last %s components in the right argument."%(axis_offset,axis_offset)
4680           DependendSymbol.__init__(self,dim=commonDim(arg0,arg1),shape=sh0+sh1,args=[arg0,arg1,axis_offset])
4681    
4682       def getMyCode(self,argstrs,format="escript"):
4683          """
4684          returns a program code that can be used to evaluate the symbol.
4685    
4686          @param argstrs: gives for each argument a string representing the argument for the evaluation.
4687          @type argstrs: C{list} of length 2 of C{str}.
4688          @param format: specifies the format to be used. At the moment only "escript", "str" and "text" are supported.
4689          @type format: C{str}
4690          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
4691          @rtype: C{str}
4692          @raise NotImplementedError: if the requested format is not available
4693          """
4694          if format=="escript" or format=="str" or format=="text":
4695             return "generalTensorTransposedProduct(%s,%s,axis_offset=%s)"%(argstrs[0],argstrs[1],argstrs[2])
4696          else:
4697             raise NotImplementedError,"%s does not provide program code for format %s."%(str(self),format)
4698    
4699       def substitute(self,argvals):
4700          """
4701          assigns new values to symbols in the definition of the symbol.
4702          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
4703    
4704          @param argvals: new values assigned to symbols
4705          @type argvals: C{dict} with keywords of type L{Symbol}.
4706          @return: result of the substitution process. Operations are executed as much as possible.
4707          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
4708          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
4709          """
4710          if argvals.has_key(self):
4711             arg=argvals[self]
4712             if self.isAppropriateValue(arg):
4713                return arg
4714             else:
4715                raise TypeError,"%s: new value is not appropriate."%str(self)
4716          else:
4717             args=self.getSubstitutedArguments(argvals)
4718             return generalTensorTransposedProduct(args[0],args[1],args[2])
4719    
4720    def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset): # this should be escript._generalTensorTransposedProduct
4721        "arg0 and arg1 are both Data objects but not neccesrily on the same function space. they could be identical!!!"
4722        return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
4723    
4724  #=========================================================  #=========================================================
4725  #  functions dealing with spatial dependency  #  functions dealing with spatial dependency

Legend:
Removed from v.775  
changed lines
  Added in v.881

  ViewVC Help
Powered by ViewVC 1.1.26