/[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 443 by gross, Fri Jan 20 06:22:38 2006 UTC revision 587 by gross, Fri Mar 10 02:26:50 2006 UTC
# Line 44  import os Line 44  import os
44  # def matchType(arg0=0.,arg1=0.):  # def matchType(arg0=0.,arg1=0.):
45  # def matchShape(arg0,arg1):  # def matchShape(arg0,arg1):
46    
 # def transpose(arg,axis=None):  
47  # def reorderComponents(arg,index):  # def reorderComponents(arg,index):
48    
49  #  #
# Line 134  def identity(shape=()): Line 133  def identity(shape=()):
133     @raise ValueError: if len(shape)>2.     @raise ValueError: if len(shape)>2.
134     """     """
135     if len(shape)>0:     if len(shape)>0:
136        out=numarray.zeros(shape+shape,numarray.Float)        out=numarray.zeros(shape+shape,numarray.Float64)
137        if len(shape)==1:        if len(shape)==1:
138            for i0 in range(shape[0]):            for i0 in range(shape[0]):
139               out[i0,i0]=1.               out[i0,i0]=1.
# Line 390  def matchType(arg0=0.,arg1=0.): Line 389  def matchType(arg0=0.,arg1=0.):
389         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
390            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
391         elif isinstance(arg1,float):         elif isinstance(arg1,float):
392            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
393         elif isinstance(arg1,int):         elif isinstance(arg1,int):
394            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
395         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
396            pass            pass
397         else:         else:
# Line 416  def matchType(arg0=0.,arg1=0.): Line 415  def matchType(arg0=0.,arg1=0.):
415         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
416            pass            pass
417         elif isinstance(arg1,float):         elif isinstance(arg1,float):
418            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
419         elif isinstance(arg1,int):         elif isinstance(arg1,int):
420            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
421         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
422            pass            pass
423         else:         else:
424            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
425      elif isinstance(arg0,float):      elif isinstance(arg0,float):
426         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
427            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
428         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
429            arg0=escript.Data(arg0,arg1.getFunctionSpace())            arg0=escript.Data(arg0,arg1.getFunctionSpace())
430         elif isinstance(arg1,float):         elif isinstance(arg1,float):
431            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
432            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
433         elif isinstance(arg1,int):         elif isinstance(arg1,int):
434            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
435            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
436         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
437            arg0=numarray.array(arg0)            arg0=numarray.array(arg0,type=numarray.Float64)
438         else:         else:
439            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
440      elif isinstance(arg0,int):      elif isinstance(arg0,int):
441         if isinstance(arg1,numarray.NumArray):         if isinstance(arg1,numarray.NumArray):
442            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
443         elif isinstance(arg1,escript.Data):         elif isinstance(arg1,escript.Data):
444            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())            arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
445         elif isinstance(arg1,float):         elif isinstance(arg1,float):
446            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
447            arg1=numarray.array(arg1)            arg1=numarray.array(arg1,type=numarray.Float64)
448         elif isinstance(arg1,int):         elif isinstance(arg1,int):
449            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
450            arg1=numarray.array(float(arg1))            arg1=numarray.array(float(arg1),type=numarray.Float64)
451         elif isinstance(arg1,Symbol):         elif isinstance(arg1,Symbol):
452            arg0=numarray.array(float(arg0))            arg0=numarray.array(float(arg0),type=numarray.Float64)
453         else:         else:
454            raise TypeError,"function: Unknown type of second argument."                raise TypeError,"function: Unknown type of second argument."    
455      else:      else:
# Line 473  def matchShape(arg0,arg1): Line 472  def matchShape(arg0,arg1):
472      sh0=pokeShape(arg0)      sh0=pokeShape(arg0)
473      sh1=pokeShape(arg1)      sh1=pokeShape(arg1)
474      if len(sh0)<len(sh):      if len(sh0)<len(sh):
475         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float)),arg1         return outer(arg0,numarray.ones(sh[len(sh0):],numarray.Float64)),arg1
476      elif len(sh1)<len(sh):      elif len(sh1)<len(sh):
477         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float))         return arg0,outer(arg1,numarray.ones(sh[len(sh1):],numarray.Float64))
478      else:      else:
479         return arg0,arg1         return arg0,arg1
480  #=========================================================  #=========================================================
# Line 601  class Symbol(object): Line 600  class Symbol(object):
600            else:            else:
601                s=pokeShape(s)+arg.getShape()                s=pokeShape(s)+arg.getShape()
602                if len(s)>0:                if len(s)>0:
603                   out.append(numarray.zeros(s),numarray.Float)                   out.append(numarray.zeros(s),numarray.Float64)
604                else:                else:
605                   out.append(a)                   out.append(a)
606         return out         return out
# Line 691  class Symbol(object): Line 690  class Symbol(object):
690         else:         else:
691            s=self.getShape()+arg.getShape()            s=self.getShape()+arg.getShape()
692            if len(s)>0:            if len(s)>0:
693               return numarray.zeros(s,numarray.Float)               return numarray.zeros(s,numarray.Float64)
694            else:            else:
695               return 0.               return 0.
696    
# Line 829  class Symbol(object): Line 828  class Symbol(object):
828         """         """
829         return power(other,self)         return power(other,self)
830    
831       def __getitem__(self,index):
832           """
833           returns the slice defined by index
834    
835           @param index: defines a
836           @type index: C{slice} or C{int} or a C{tuple} of them
837           @return: a S{Symbol} representing the slice defined by index
838           @rtype: L{DependendSymbol}
839           """
840           return GetSlice_Symbol(self,index)
841    
842  class DependendSymbol(Symbol):  class DependendSymbol(Symbol):
843     """     """
844     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.     DependendSymbol extents L{Symbol} by modifying the == operator to allow two instances to be equal.
# Line 879  class DependendSymbol(Symbol): Line 889  class DependendSymbol(Symbol):
889  #=========================================================  #=========================================================
890  #  Unary operations prserving the shape  #  Unary operations prserving the shape
891  #========================================================  #========================================================
892    class GetSlice_Symbol(DependendSymbol):
893       """
894       L{Symbol} representing getting a slice for a L{Symbol}
895       """
896       def __init__(self,arg,index):
897          """
898          initialization of wherePositive L{Symbol} with argument arg
899          @param arg: argument
900          @type arg: L{Symbol}.
901          @param index: defines index
902          @type index: C{slice} or C{int} or a C{tuple} of them
903          @raises IndexError: if length of index is larger than rank of arg or a index start or stop is out of range
904          @raises ValueError: if a step is given
905          """
906          if not isinstance(index,tuple): index=(index,)
907          if len(index)>arg.getRank():
908               raise IndexError,"GetSlice_Symbol: index out of range."
909          sh=()
910          index2=()
911          for i in range(len(index)):
912             ix=index[i]
913             if isinstance(ix,int):
914                if ix<0 or ix>=arg.getShape()[i]:
915                   raise ValueError,"GetSlice_Symbol: index out of range."
916                index2=index2+(ix,)
917             else:
918               if not ix.step==None:
919                 raise ValueError,"GetSlice_Symbol: steping is not supported."
920               if ix.start==None:
921                  s=0
922               else:
923                  s=ix.start
924               if ix.stop==None:
925                  e=arg.getShape()[i]
926               else:
927                  e=ix.stop
928                  if e>arg.getShape()[i]:
929                     raise IndexError,"GetSlice_Symbol: index out of range."
930               index2=index2+(slice(s,e),)
931               if e>s:
932                   sh=sh+(e-s,)
933               elif s>e:
934                   raise IndexError,"GetSlice_Symbol: slice start must be less or equal slice end"
935          for i in range(len(index),arg.getRank()):
936              index2=index2+(slice(0,arg.getShape()[i]),)
937              sh=sh+(arg.getShape()[i],)
938          super(GetSlice_Symbol, self).__init__(args=[arg,index2],shape=sh,dim=arg.getDim())
939    
940       def getMyCode(self,argstrs,format="escript"):
941          """
942          returns a program code that can be used to evaluate the symbol.
943    
944          @param argstrs: gives for each argument a string representing the argument for the evaluation.
945          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
946          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
947          @type format: C{str}
948          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
949          @rtype: C{str}
950          @raise: NotImplementedError: if the requested format is not available
951          """
952          if format=="escript" or format=="str"  or format=="text":
953             return "%s.__getitem__(%s)"%(argstrs[0],argstrs[1])
954          else:
955             raise NotImplementedError,"GetItem_Symbol does not provide program code for format %s."%format
956    
957       def substitute(self,argvals):
958          """
959          assigns new values to symbols in the definition of the symbol.
960          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
961    
962          @param argvals: new values assigned to symbols
963          @type argvals: C{dict} with keywords of type L{Symbol}.
964          @return: result of the substitution process. Operations are executed as much as possible.
965          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
966          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
967          """
968          if argvals.has_key(self):
969             arg=argvals[self]
970             if self.isAppropriateValue(arg):
971                return arg
972             else:
973                raise TypeError,"%s: new value is not appropriate."%str(self)
974          else:
975             args=self.getSubstitutedArguments(argvals)
976             arg=args[0]
977             index=args[1]
978             return arg.__getitem__(index)
979    
980  def log10(arg):  def log10(arg):
981     """     """
982     returns base-10 logarithm of argument arg     returns base-10 logarithm of argument arg
# Line 911  def wherePositive(arg): Line 1009  def wherePositive(arg):
1009     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1010     """     """
1011     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1012        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1013        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1014        return out        return out
1015     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1016        return arg._wherePositive()        return arg._wherePositive()
# Line 993  def whereNegative(arg): Line 1091  def whereNegative(arg):
1091     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1092     """     """
1093     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1094        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1095        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1096        return out        return out
1097     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1098        return arg._whereNegative()        return arg._whereNegative()
# Line 1075  def whereNonNegative(arg): Line 1173  def whereNonNegative(arg):
1173     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1174     """     """
1175     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1176        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1177        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1178        return out        return out
1179     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1180        return arg._whereNonNegative()        return arg._whereNonNegative()
# Line 1105  def whereNonPositive(arg): Line 1203  def whereNonPositive(arg):
1203     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1204     """     """
1205     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1206        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(arg,numarray.zeros(arg.shape,numarray.Float64))*1.
1207        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1208        return out        return out
1209     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1210        return arg._whereNonPositive()        return arg._whereNonPositive()
# Line 1137  def whereZero(arg,tol=0.): Line 1235  def whereZero(arg,tol=0.):
1235     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1236     """     """
1237     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1238        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.less_equal(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1239        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1240        return out        return out
1241     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1242        if tol>0.:        if tol>0.:
# Line 1220  def whereNonZero(arg,tol=0.): Line 1318  def whereNonZero(arg,tol=0.):
1318     @raises TypeError: if the type of the argument is not expected.     @raises TypeError: if the type of the argument is not expected.
1319     """     """
1320     if isinstance(arg,numarray.NumArray):     if isinstance(arg,numarray.NumArray):
1321        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float))*1.        out=numarray.greater(abs(arg)-tol,numarray.zeros(arg.shape,numarray.Float64))*1.
1322        if isinstance(out,float): out=numarray.array(out)        if isinstance(out,float): out=numarray.array(out,type=numarray.Float64)
1323        return out        return out
1324     elif isinstance(arg,escript.Data):     elif isinstance(arg,escript.Data):
1325        if tol>0.:        if tol>0.:
# Line 2882  def trace(arg,axis_offset=0): Line 2980  def trace(arg,axis_offset=0):
2980        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
2981          raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"trace: dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)
2982        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))
2983        out=numarray.zeros([s1,s2],numarray.Float)        out=numarray.zeros([s1,s2],numarray.Float64)
2984        for i1 in range(s1):        for i1 in range(s1):
2985          for i2 in range(s2):          for i2 in range(s2):
2986              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]              for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
# Line 3012  class Trace_Symbol(DependendSymbol): Line 3110  class Trace_Symbol(DependendSymbol):
3110        else:        else:
3111           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])           return trace(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3112    
3113    def transpose(arg,axis_offset=None):
3114       """
3115       returns the transpose of arg by swaping the first axis_offset and the last rank-axis_offset components.
3116    
3117       @param arg: argument
3118       @type arg: L{escript.Data}, L{Symbol}, L{numarray.NumArray}, C{float}, C{int}
3119       @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.
3120                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3121       @type axis_offset: C{int}
3122       @return: transpose of arg
3123       @rtype: L{escript.Data}, L{Symbol}, L{numarray.NumArray},C{float}, C{int} depending on the type of arg.
3124       """
3125       if isinstance(arg,numarray.NumArray):
3126          if axis_offset==None: axis_offset=int(arg.rank/2)
3127          return numarray.transpose(arg,axes=range(axis_offset,arg.rank)+range(0,axis_offset))
3128       elif isinstance(arg,escript.Data):
3129          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3130          return escript_transpose(arg,axis_offset)
3131       elif isinstance(arg,float):
3132          if not ( axis_offset==0 or axis_offset==None):
3133            raise ValueError,"transpose: axis_offset must be 0 for float argument"
3134          return arg
3135       elif isinstance(arg,int):
3136          if not ( axis_offset==0 or axis_offset==None):
3137            raise ValueError,"transpose: axis_offset must be 0 for int argument"
3138          return float(arg)
3139       elif isinstance(arg,Symbol):
3140          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3141          return Transpose_Symbol(arg,axis_offset)
3142       else:
3143          raise TypeError,"transpose: Unknown argument type."
3144    
3145    def escript_transpose(arg,axis_offset): # this should be escript._transpose
3146          "arg si a Data objects!!!"
3147          r=arg.getRank()
3148          if axis_offset<0 or axis_offset>r:
3149            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3150          s=arg.getShape()
3151          s_out=s[axis_offset:]+s[:axis_offset]
3152          out=escript.Data(0.,s_out,arg.getFunctionSpace())
3153          if r==4:
3154             if axis_offset==1:
3155                for i0 in range(s_out[0]):
3156                   for i1 in range(s_out[1]):
3157                      for i2 in range(s_out[2]):
3158                         for i3 in range(s_out[3]):
3159                             out[i0,i1,i2,i3]=arg[i3,i0,i1,i2]
3160             elif axis_offset==2:
3161                for i0 in range(s_out[0]):
3162                   for i1 in range(s_out[1]):
3163                      for i2 in range(s_out[2]):
3164                         for i3 in range(s_out[3]):
3165                             out[i0,i1,i2,i3]=arg[i2,i3,i0,i1]
3166             elif axis_offset==3:
3167                for i0 in range(s_out[0]):
3168                   for i1 in range(s_out[1]):
3169                      for i2 in range(s_out[2]):
3170                         for i3 in range(s_out[3]):
3171                             out[i0,i1,i2,i3]=arg[i1,i2,i3,i0]
3172             else:
3173                for i0 in range(s_out[0]):
3174                   for i1 in range(s_out[1]):
3175                      for i2 in range(s_out[2]):
3176                         for i3 in range(s_out[3]):
3177                             out[i0,i1,i2,i3]=arg[i0,i1,i2,i3]
3178          elif r==3:
3179             if axis_offset==1:
3180                for i0 in range(s_out[0]):
3181                   for i1 in range(s_out[1]):
3182                      for i2 in range(s_out[2]):
3183                             out[i0,i1,i2]=arg[i2,i0,i1]
3184             elif axis_offset==2:
3185                for i0 in range(s_out[0]):
3186                   for i1 in range(s_out[1]):
3187                      for i2 in range(s_out[2]):
3188                             out[i0,i1,i2]=arg[i1,i2,i0]
3189             else:
3190                for i0 in range(s_out[0]):
3191                   for i1 in range(s_out[1]):
3192                      for i2 in range(s_out[2]):
3193                             out[i0,i1,i2]=arg[i0,i1,i2]
3194          elif r==2:
3195             if axis_offset==1:
3196                for i0 in range(s_out[0]):
3197                   for i1 in range(s_out[1]):
3198                             out[i0,i1]=arg[i1,i0]
3199             else:
3200                for i0 in range(s_out[0]):
3201                   for i1 in range(s_out[1]):
3202                             out[i0,i1]=arg[i0,i1]
3203          elif r==1:
3204              for i0 in range(s_out[0]):
3205                   out[i0]=arg[i0]
3206          elif r==0:
3207                 out=arg+0.
3208          return out
3209    class Transpose_Symbol(DependendSymbol):
3210       """
3211       L{Symbol} representing the result of the transpose function
3212       """
3213       def __init__(self,arg,axis_offset=None):
3214          """
3215          initialization of transpose L{Symbol} with argument arg
3216    
3217          @param arg: argument of function
3218          @type arg: L{Symbol}.
3219           @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.
3220                           if axis_offset is not present C{int(r/2)} where r is the rank of arg is used.
3221          @type axis_offset: C{int}
3222          """
3223          if axis_offset==None: axis_offset=int(arg.getRank()/2)
3224          if axis_offset<0 or axis_offset>arg.getRank():
3225            raise ValueError,"escript_transpose: axis_offset must be between 0 and %s"%r
3226          s=arg.getShape()
3227          super(Transpose_Symbol,self).__init__(args=[arg,axis_offset],shape=s[axis_offset:]+s[:axis_offset],dim=arg.getDim())
3228    
3229       def getMyCode(self,argstrs,format="escript"):
3230          """
3231          returns a program code that can be used to evaluate the symbol.
3232    
3233          @param argstrs: gives for each argument a string representing the argument for the evaluation.
3234          @type argstrs: C{str} or a C{list} of length 1 of C{str}.
3235          @param format: specifies the format to be used. At the moment only "escript" ,"text" and "str" are supported.
3236          @type format: C{str}
3237          @return: a piece of program code which can be used to evaluate the expression assuming the values for the arguments are available.
3238          @rtype: C{str}
3239          @raise: NotImplementedError: if the requested format is not available
3240          """
3241          if format=="escript" or format=="str"  or format=="text":
3242             return "transpose(%s,axis_offset=%s)"%(argstrs[0],argstrs[1])
3243          else:
3244             raise NotImplementedError,"Transpose_Symbol does not provide program code for format %s."%format
3245    
3246       def substitute(self,argvals):
3247          """
3248          assigns new values to symbols in the definition of the symbol.
3249          The method replaces the L{Symbol} u by argvals[u] in the expression defining this object.
3250    
3251          @param argvals: new values assigned to symbols
3252          @type argvals: C{dict} with keywords of type L{Symbol}.
3253          @return: result of the substitution process. Operations are executed as much as possible.
3254          @rtype: L{escript.Symbol}, C{float}, L{escript.Data}, L{numarray.NumArray} depending on the degree of substitution
3255          @raise TypeError: if a value for a L{Symbol} cannot be substituted.
3256          """
3257          if argvals.has_key(self):
3258             arg=argvals[self]
3259             if self.isAppropriateValue(arg):
3260                return arg
3261             else:
3262                raise TypeError,"%s: new value is not appropriate."%str(self)
3263          else:
3264             arg=self.getSubstitutedArguments(argvals)
3265             return transpose(arg[0],axis_offset=arg[1])
3266    
3267       def diff(self,arg):
3268          """
3269          differential of this object
3270    
3271          @param arg: the derivative is calculated with respect to arg
3272          @type arg: L{escript.Symbol}
3273          @return: derivative with respect to C{arg}
3274          @rtype: typically L{Symbol} but other types such as C{float}, L{escript.Data}, L{numarray.NumArray}  are possible.
3275          """
3276          if arg==self:
3277             return identity(self.getShape())
3278          else:
3279             return transpose(self.getDifferentiatedArguments(arg)[0],axis_offset=self.getArgument()[1])
3280    def symmetric(arg):
3281        """
3282        returns the symmetric part of the square matrix arg. This is (arg+transpose(arg))/2
3283    
3284        @param arg: square matrix. Must have rank 2 or 4 and be square.
3285        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3286        @return: symmetric part of arg
3287        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3288        """
3289        if isinstance(arg,numarray.NumArray):
3290          if arg.rank==2:
3291            if not (arg.shape[0]==arg.shape[1]):
3292               raise ValueError,"symmetric: argument must be square."
3293          elif arg.rank==4:
3294            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3295               raise ValueError,"symmetric: argument must be square."
3296          else:
3297            raise ValueError,"symmetric: rank 2 or 4 is required."
3298          return (arg+transpose(arg))/2
3299        elif isinstance(arg,escript.Data):
3300          return escript_symmetric(arg)
3301        elif isinstance(arg,float):
3302          return arg
3303        elif isinstance(arg,int):
3304          return float(arg)
3305        elif isinstance(arg,Symbol):
3306          if arg.getRank()==2:
3307            if not (arg.getShape()[0]==arg.getShape()[1]):
3308               raise ValueError,"symmetric: argument must be square."
3309          elif arg.getRank()==4:
3310            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3311               raise ValueError,"symmetric: argument must be square."
3312          else:
3313            raise ValueError,"symmetric: rank 2 or 4 is required."
3314          return (arg+transpose(arg))/2
3315        else:
3316          raise TypeError,"symmetric: Unknown argument type."
3317    
3318    def escript_symmetric(arg): # this should be implemented in c++
3319          if arg.getRank()==2:
3320            if not (arg.getShape()[0]==arg.getShape()[1]):
3321               raise ValueError,"escript_symmetric: argument must be square."
3322            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3323            for i0 in range(arg.getShape()[0]):
3324               for i1 in range(arg.getShape()[1]):
3325                  out[i0,i1]=(arg[i0,i1]+arg[i1,i0])/2.
3326          elif arg.getRank()==4:
3327            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3328               raise ValueError,"escript_symmetric: argument must be square."
3329            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3330            for i0 in range(arg.getShape()[0]):
3331               for i1 in range(arg.getShape()[1]):
3332                  for i2 in range(arg.getShape()[2]):
3333                     for i3 in range(arg.getShape()[3]):
3334                         out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]+arg[i2,i3,i0,i1])/2.
3335          else:
3336            raise ValueError,"escript_symmetric: rank 2 or 4 is required."
3337          return out
3338    
3339    def nonsymmetric(arg):
3340        """
3341        returns the nonsymmetric part of the square matrix arg. This is (arg-transpose(arg))/2
3342    
3343        @param arg: square matrix. Must have rank 2 or 4 and be square.
3344        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3345        @return: nonsymmetric part of arg
3346        @rtype: L{numarray.NumArray}, L{escript.Data}, L{Symbol} depending on the input
3347        """
3348        if isinstance(arg,numarray.NumArray):
3349          if arg.rank==2:
3350            if not (arg.shape[0]==arg.shape[1]):
3351               raise ValueError,"nonsymmetric: argument must be square."
3352          elif arg.rank==4:
3353            if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
3354               raise ValueError,"nonsymmetric: argument must be square."
3355          else:
3356            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3357          return (arg-transpose(arg))/2
3358        elif isinstance(arg,escript.Data):
3359          return escript_nonsymmetric(arg)
3360        elif isinstance(arg,float):
3361          return arg
3362        elif isinstance(arg,int):
3363          return float(arg)
3364        elif isinstance(arg,Symbol):
3365          if arg.getRank()==2:
3366            if not (arg.getShape()[0]==arg.getShape()[1]):
3367               raise ValueError,"nonsymmetric: argument must be square."
3368          elif arg.getRank()==4:
3369            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3370               raise ValueError,"nonsymmetric: argument must be square."
3371          else:
3372            raise ValueError,"nonsymmetric: rank 2 or 4 is required."
3373          return (arg-transpose(arg))/2
3374        else:
3375          raise TypeError,"nonsymmetric: Unknown argument type."
3376    
3377    def escript_nonsymmetric(arg): # this should be implemented in c++
3378          if arg.getRank()==2:
3379            if not (arg.getShape()[0]==arg.getShape()[1]):
3380               raise ValueError,"escript_nonsymmetric: argument must be square."
3381            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3382            for i0 in range(arg.getShape()[0]):
3383               for i1 in range(arg.getShape()[1]):
3384                  out[i0,i1]=(arg[i0,i1]-arg[i1,i0])/2.
3385          elif arg.getRank()==4:
3386            if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
3387               raise ValueError,"escript_nonsymmetric: argument must be square."
3388            out=escript.Data(0.,arg.getShape(),arg.getFunctionSpace())
3389            for i0 in range(arg.getShape()[0]):
3390               for i1 in range(arg.getShape()[1]):
3391                  for i2 in range(arg.getShape()[2]):
3392                     for i3 in range(arg.getShape()[3]):
3393                         out[i0,i1,i2,i3]=(arg[i0,i1,i2,i3]-arg[i2,i3,i0,i1])/2.
3394          else:
3395            raise ValueError,"escript_nonsymmetric: rank 2 or 4 is required."
3396          return out
3397    
3398    
3399  def inverse(arg):  def inverse(arg):
3400      """      """
3401      returns the inverse of the square matrix arg.      returns the inverse of the square matrix arg.
3402    
3403      @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.
3404      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}      @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3405      @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 matrixmul(inverse(arg),arg) almost equal to kronecker(arg.getShape()[0])
3406      @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
3407        @remark: for L{escript.Data} objects the dimension is restricted to 3.
3408      """      """
3409      if isinstance(arg,numarray.NumArray):      if isinstance(arg,numarray.NumArray):
3410        return numarray.linear_algebra.inverse(arg)        return numarray.linear_algebra.inverse(arg)
# Line 3154  class Inverse_Symbol(DependendSymbol): Line 3539  class Inverse_Symbol(DependendSymbol):
3539           return identity(self.getShape())           return identity(self.getShape())
3540        else:        else:
3541           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)           return -matrixmult(matrixmult(self,self.getDifferentiatedArguments(arg)[0]),self)
3542    
3543    def eigenvalues(arg):
3544        """
3545        returns the eigenvalues of the square matrix arg.
3546    
3547        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3548                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3549        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3550        @return: the eigenvalues in increasing order.
3551        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3552        @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3553        """
3554        if isinstance(arg,numarray.NumArray):
3555          out=numarray.linear_algebra.eigenvalues((arg+numarray.transpose(arg))/2.)
3556          out.sort()
3557          return out
3558        elif isinstance(arg,escript.Data):
3559          return arg._eigenvalues()
3560        elif isinstance(arg,Symbol):
3561          if not arg.getRank()==2:
3562            raise ValueError,"eigenvalues: argument must have rank 2"
3563          s=arg.getShape()      
3564          if not s[0] == s[1]:
3565            raise ValueError,"eigenvalues: argument must be a square matrix."
3566          if s[0]==1:
3567              return arg[0]
3568          elif s[0]==2:
3569              arg1=symmetric(arg)
3570              A11=arg1[0,0]
3571              A12=arg1[0,1]
3572              A22=arg1[1,1]
3573              trA=(A11+A22)/2.
3574              A11-=trA
3575              A22-=trA
3576              s=sqrt(A12**2-A11*A22)
3577              return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3578          elif s[0]==3:
3579              arg1=symmetric(arg)
3580              A11=arg1[0,0]
3581              A12=arg1[0,1]
3582              A22=arg1[1,1]
3583              A13=arg1[0,2]
3584              A23=arg1[1,2]
3585              A33=arg1[2,2]
3586              trA=(A11+A22+A33)/3.
3587              A11-=trA
3588              A22-=trA
3589              A33-=trA
3590              A13_2=A13**2
3591              A23_2=A23**2
3592              A12_2=A12**2
3593              p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3594              q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3595              sq_p=sqrt(p/3.)
3596              alpha_3=acos(clip(-q*(sq_p+whereZero(p,0.)*1.e-15)**(-3.)/2.,-1.,1.))/3.  # whereZero is protection against divison by zero
3597              sq_p*=2.
3598              f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3599               -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
3600               -cos(alpha_3-numarray.pi/3.)*numarray.array([1.,0.,0.],type=numarray.Float64)
3601              return trA+sq_p*f
3602          else:
3603             raise TypeError,"eigenvalues: only matrix dimensions 1,2,3 are supported right now."
3604        elif isinstance(arg,float):
3605          return arg
3606        elif isinstance(arg,int):
3607          return float(arg)
3608        else:
3609          raise TypeError,"eigenvalues: Unknown argument type."
3610    
3611    def eigenvalues_and_eigenvectors(arg):
3612        """
3613        returns the eigenvalues of the square matrix arg.
3614    
3615        @param arg: square matrix. Must have rank 2 and the first and second dimension must be equal.
3616                    arg must be symmetric, ie. transpose(arg)==arg (this is not checked).
3617        @type arg: L{numarray.NumArray}, L{escript.Data}, L{Symbol}
3618        @return: the eigenvalues in increasing order.
3619        @rtype: L{numarray.NumArray},L{escript.Data}, L{Symbol} depending on the input.
3620        @remark: for L{escript.Data} and L{Symbol} objects the dimension is restricted to 3.
3621        """
3622        if isinstance(arg,numarray.NumArray):
3623          raise TypeError,"eigenvalues_and_eigenvectors is not supporting numarray arguments"
3624        elif isinstance(arg,escript.Data):
3625          return arg._eigenvalues_and_eigenvectors()
3626        elif isinstance(arg,Symbol):
3627          raise TypeError,"eigenvalues_and_eigenvectors is not supporting Symbol arguments"
3628        elif isinstance(arg,float):
3629          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3630        elif isinstance(arg,int):
3631          return (numarray.array([[arg]],numarray.Float),numarray.ones((1,1),numarray.Float))
3632        else:
3633          raise TypeError,"eigenvalues: Unknown argument type."
3634  #=======================================================  #=======================================================
3635  #  Binary operations:  #  Binary operations:
3636  #=======================================================  #=======================================================
# Line 3272  def mult(arg0,arg1): Line 3749  def mult(arg0,arg1):
3749         """         """
3750         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3751         if testForZero(args[0]) or testForZero(args[1]):         if testForZero(args[0]) or testForZero(args[1]):
3752            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3753         else:         else:
3754            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :            if isinstance(args[0],Symbol) or isinstance(args[1],Symbol) :
3755                return Mult_Symbol(args[0],args[1])                return Mult_Symbol(args[0],args[1])
# Line 3372  def quotient(arg0,arg1): Line 3849  def quotient(arg0,arg1):
3849         """         """
3850         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3851         if testForZero(args[0]):         if testForZero(args[0]):
3852            return numarray.zeros(pokeShape(args[0]),numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3853         elif isinstance(args[0],Symbol):         elif isinstance(args[0],Symbol):
3854            if isinstance(args[1],Symbol):            if isinstance(args[1],Symbol):
3855               return Quotient_Symbol(args[0],args[1])               return Quotient_Symbol(args[0],args[1])
# Line 3478  def power(arg0,arg1): Line 3955  def power(arg0,arg1):
3955         """         """
3956         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3957         if testForZero(args[0]):         if testForZero(args[0]):
3958            return numarray.zeros(args[0],numarray.Float)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3959         elif testForZero(args[1]):         elif testForZero(args[1]):
3960            return numarray.ones(args[0],numarray.Float)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3961         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3962            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3963         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 3777  def generalTensorProduct(arg0,arg1,axis_ Line 4254  def generalTensorProduct(arg0,arg1,axis_
4254             for i in sh1[:axis_offset]: d01*=i             for i in sh1[:axis_offset]: d01*=i
4255             arg0_c.resize((d0,d01))             arg0_c.resize((d0,d01))
4256             arg1_c.resize((d01,d1))             arg1_c.resize((d01,d1))
4257             out=numarray.zeros((d0,d1),numarray.Float)             out=numarray.zeros((d0,d1),numarray.Float64)
4258             for i0 in range(d0):             for i0 in range(d0):
4259                      for i1 in range(d1):                      for i1 in range(d1):
4260                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])                           out[i0,i1]=numarray.sum(arg0_c[i0,:]*arg1_c[:,i1])
# Line 4224  def L2(arg): Line 4701  def L2(arg):
4701      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))
4702  #=============================  #=============================
4703  #  #
 # wrapper for various functions: if the argument has attribute the function name  
 # as an argument it calls the corresponding methods. Otherwise the corresponding  
 # numarray function is called.  
   
 # functions involving the underlying Domain:  
   
 def transpose(arg,axis=None):  
     """  
     Returns the transpose of the Data object arg.  
   
     @param arg:  
     """  
     if axis==None:  
        r=0  
        if hasattr(arg,"getRank"): r=arg.getRank()  
        if hasattr(arg,"rank"): r=arg.rank  
        axis=r/2  
     if isinstance(arg,Symbol):  
        return Transpose_Symbol(arg,axis=r)  
     if isinstance(arg,escript.Data):  
        # hack for transpose  
        r=arg.getRank()  
        if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"  
        s=arg.getShape()  
        out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())  
        for i in range(s[0]):  
           for j in range(s[1]):  
              out[j,i]=arg[i,j]  
        return out  
        # end hack for transpose  
        return arg.transpose(axis)  
     else:  
        return numarray.transpose(arg,axis=axis)  
   
   
4704    
4705  def reorderComponents(arg,index):  def reorderComponents(arg,index):
4706      """      """
4707      resorts the component of arg according to index      resorts the component of arg according to index
4708    
4709      """      """
4710      pass      raise NotImplementedError
4711  #  #
4712  # $Log: util.py,v $  # $Log: util.py,v $
4713  # Revision 1.14.2.16  2005/10/19 06:09:57  gross  # Revision 1.14.2.16  2005/10/19 06:09:57  gross
Legend:
Removed from v.443  
changed lines
  Added in v.587
  ViewVC Help
Powered by ViewVC 1.1.26