/[escript]/branches/symbolic_from_3470/escript/py_src/util.py
ViewVC logotype

Annotation of /branches/symbolic_from_3470/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3631 - (hide annotations)
Thu Oct 20 03:12:46 2011 UTC (7 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 104513 byte(s)
Fix compatibility with latest sympy version.

1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1312 #
12 ksteube 1809 ########################################################
13 jgs 82
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 jgs 82 """
23 jgs 149 Utility functions for escript
24 jgs 123
25 jfenwick 2625 :var __author__: name of author
26     :var __copyright__: copyrights
27     :var __license__: licence agreement
28     :var __url__: url entry point on documentation
29     :var __version__: version
30     :var __date__: date of the version
31     :var EPSILON: smallest positive value with 1.<1.+EPSILON
32     :var DBLE_MAX: largest positive float
33 jgs 82 """
34 ksteube 1809
35 gross 290 __author__="Lutz Gross, l.gross@uq.edu.au"
36 jgs 82
37 gross 290
38     import math
39 jfenwick 2455 import numpy
40 jgs 102 import escript
41 jgs 150 import os
42 ksteube 813 from esys.escript import C_GeneralTensorProduct
43 gross 2415 from esys.escript import getVersion, getMPIRankWorld, getMPIWorldMax
44 ksteube 1561 from esys.escript import printParallelThreadCounts
45 jfenwick 2142 from esys.escript import listEscriptParams
46 jfenwick 3059 from esys.escript.escriptcpp import Data, _saveDataCSV, _condEval
47 caltinay 3631 from esys.escript.symbolic import *
48 jgs 124
49 jgs 123 #=========================================================
50 gross 290 # some helpers:
51 jgs 123 #=========================================================
52 ksteube 1312 def getEpsilon():
53 gross 2100 return escript.getMachinePrecision()
54 ksteube 1312 EPSILON=getEpsilon()
55    
56 gross 2100 def getMaxFloat():
57     return escript.getMaxFloat()
58     DBLE_MAX=getMaxFloat()
59    
60 gross 1042 def getTagNames(domain):
61     """
62 caltinay 2169 Returns a list of tag names used by the domain.
63 gross 1042
64 jfenwick 2625 :param domain: a domain object
65     :type domain: `escript.Domain`
66     :return: a list of tag names used by the domain
67     :rtype: ``list`` of ``str``
68 gross 1042 """
69     return [n.strip() for n in domain.showTagNames().split(",") ]
70 gross 1044
71     def insertTagNames(domain,**kwargs):
72     """
73 caltinay 2169 Inserts tag names into the domain.
74 gross 1044
75 jfenwick 2625 :param domain: a domain object
76     :type domain: ``escript.Domain``
77     :keyword <tag_name>: tag key assigned to <tag_name>
78     :type <tag_name>: ``int``
79 gross 1044 """
80     for k in kwargs:
81     domain.setTagMap(k,kwargs[k])
82    
83     def insertTaggedValues(target,**kwargs):
84     """
85 caltinay 2169 Inserts tagged values into the target using tag names.
86 gross 1044
87 jfenwick 2625 :param target: data to be filled by tagged values
88     :type target: `escript.Data`
89     :keyword <tag_name>: value to be used for <tag_name>
90     :type <tag_name>: ``float`` or ``numpy.ndarray``
91     :return: ``target``
92     :rtype: `escript.Data`
93 gross 1044 """
94     for k in kwargs:
95     target.setTaggedValue(k,kwargs[k])
96     return target
97    
98 jfenwick 3368
99     def interpolateTable(tab, dat, start, step, undef=1.e50, check_boundaries=False):
100     try:
101     dim=len(start)
102     except TypeError:
103     start=(start,)
104     dim=1
105     try:
106     slen=len(step)
107     except TypeError:
108     step=(step,)
109     slen=1
110     if dim<1 or dim>3:
111     raise ValueError("Length of start list must be between 1 and 3.")
112     if dim!=slen:
113     raise ValueError("Length of start and step must be the same.")
114     dshape=dat.getShape()
115     if len(dshape)==0:
116     datdim=0
117     firstdim=dat
118     else:
119     datdim=dshape[0]
120     firstdim=dat[0]
121     #So now we know firstdim is a scalar
122     if (dim==1 and datdim>1) or (dim>1 and datdim!=dim):
123     print dim, datdim
124     raise ValueError("The dimension of dat must be equal to the length of start.")
125     if dim==3:
126     d1=dat[1]
127     d2=dat[2]
128     return firstdim._interpolateTable3d(tab, start[0], step[0], d1, start[1], step[1], d2, start[2], step[2], undef, check_boundaries)
129     if dim==2:
130     d1=dat[1]
131     return d1.interpolateTable(tab, start[0], step[0], firstdim, start[1], step[1], undef, check_boundaries)
132     else:
133     return firstdim.interpolateTable(tab, start[0], step[0], undef, check_boundaries)
134    
135    
136 jfenwick 2635 def saveDataCSV(filename, append=False, sep=", ", csep="_", **data):
137     """
138     Writes `Data` objects to a csv file.
139     These objects must have compatible FunctionSpaces. ie it must be possible to
140     interpolate all data to one `FunctionSpace`.
141 jfenwick 3360
142     :param filename: file to save data to.
143     :type filename: ``string``
144     :param append: If ``True``, then open file at end rather than beginning
145     :type append: ``bool``
146     :param sep: separator between fields
147     :type sep: ``string``
148     :param csep: separator for components of rank2 and above eg ('_' -> c0_1)
149    
150     The keyword args are Data objects to save.
151     If a scalar `Data` object is passed with the name ``mask``, then only
152     samples which correspond to positive values in ``mask`` will be output.
153 jfenwick 2635 Example::
154 jfenwick 3360
155     s=Scalar(..)
156     v=Vector(..)
157     t=Tensor(..)
158 gross 2863 f=float()
159 jfenwick 3360 saveDataCSV("f.csv",a=s, b=v, c=t, d=f)
160 jfenwick 2635
161     Will result in a file
162    
163 gross 2863 a, b0, b1, c0_0, c0_1, .., c1_0 d
164     1.0, 1.5, 2.7, 3.1, 3.4, .., 0.89 0.0
165     0.9, 8.7, 1.9, 3.4, 7.8, .., 1.21 0.0
166 jfenwick 2635
167     The first line is a header, the remaining lines give the values.
168 jfenwick 3360
169 jfenwick 2635 """
170 gross 2863 # find a function space:
171     fs = None
172     for n,d in data.items():
173     if isinstance(d, Data): fs=d.getFunctionSpace()
174     if fs == None:
175     raise ValueError, "saveDataCSV: there must be at least one Data object in the argument list."
176    
177 jfenwick 2635 new_data={}
178     for n,d in data.items():
179 gross 2863 if isinstance(d, Data):
180     new_data[n]=d
181     else:
182     try:
183     new_data[n]=Data(d,fs)
184     except:
185     raise ValueError, "saveDataCSV: unknown non-data argument type for %s"%(str(n))
186 jfenwick 2635 _saveDataCSV(filename, new_data,sep, csep, append)
187    
188 gross 2421 def saveVTK(filename,domain=None, metadata=None, metadata_schema=None, **data):
189 jgs 150 """
190 caltinay 3344 Deprecated. See esys.weipa.saveVTK().
191     """
192 jgs 123
193 caltinay 3344 import warnings
194     msg = "esys.escript.util.saveVTK is deprecated. Use the esys.weipa module."
195     warnings.warn(msg, DeprecationWarning, stacklevel=2)
196 jgs 123
197 caltinay 3344 from esys import weipa
198     weipa.saveVTK(filename, domain, metadata, metadata_schema, **data)
199 jgs 150
200 jgs 153 def saveDX(filename,domain=None,**data):
201 jgs 149 """
202 jfenwick 2625 Writes `Data` objects into a file using the OpenDX file format.
203 jgs 149
204 gross 720 Example::
205 jgs 150
206 caltinay 2169 tmp=Scalar(..)
207     v=Vector(..)
208     saveDX("solution.dx", temperature=tmp, velocity=v)
209 jgs 150
210 jfenwick 2625 ``tmp`` and ``v`` are written into "solution.dx" where ``tmp`` is named
211     "temperature" and ``v`` is named "velocity".
212 jgs 153
213 jfenwick 2625 :param filename: file name of the output file
214     :type filename: ``str``
215     :param domain: domain of the `Data` objects. If not specified, the domain
216     of the given `Data` objects is used.
217     :type domain: `escript.Domain`
218     :keyword <name>: writes the assigned value to the DX file using <name> as
219 caltinay 2169 identifier. The identifier can be used to select the data
220     set when data are imported into DX.
221 jfenwick 2625 :type <name>: `Data` object
222     :note: The data objects have to be defined on the same domain. They may not
223     be in the same `FunctionSpace` but one cannot expect that all
224     `FunctionSpace` s can be mixed. Typically, data on the boundary and
225 caltinay 2169 data on the interior cannot be mixed.
226 jgs 149 """
227 ksteube 1312 new_data={}
228     for n,d in data.items():
229 caltinay 2169 if not d.isEmpty():
230     fs=d.getFunctionSpace()
231 ksteube 1312 domain2=fs.getDomain()
232     if fs == escript.Solution(domain2):
233     new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
234     elif fs == escript.ReducedSolution(domain2):
235     new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
236     elif fs == escript.ContinuousFunction(domain2):
237     new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
238     else:
239     new_data[n]=d
240     if domain==None: domain=domain2
241 jgs 153 if domain==None:
242 caltinay 2169 raise ValueError,"saveDX: no domain detected."
243 ksteube 1312 domain.saveDX(filename,new_data)
244 jgs 154
245 caltinay 2599 def saveESD(datasetName, dataDir=".", domain=None, timeStep=0, deltaT=1, dynamicMesh=0, **data):
246 caltinay 2188 """
247 jfenwick 2625 Saves `Data` objects to files and creates an I{escript dataset} (ESD) file
248 caltinay 2188 for convenient processing/visualisation.
249    
250 caltinay 2526 Single timestep example::
251 caltinay 2188
252     tmp = Scalar(..)
253     v = Vector(..)
254     saveESD("solution", "data", temperature=tmp, velocity=v)
255    
256 caltinay 2526 Time series example::
257    
258     while t < t_end:
259     tmp = Scalar(..)
260     v = Vector(..)
261     # save every 10 timesteps
262     if t % 10 == 0:
263     saveESD("solution", "data", timeStep=t, deltaT=10, temperature=tmp, velocity=v)
264     t = t + 1
265    
266 caltinay 2188 tmp, v and the domain are saved in native format in the "data"
267     directory and the file "solution.esd" is created that refers to tmp by
268     the name "temperature" and to v by the name "velocity".
269    
270 jfenwick 2625 :param datasetName: name of the dataset, used to name the ESD file
271     :type datasetName: ``str``
272     :param dataDir: optional directory where the data files should be saved
273     :type dataDir: ``str``
274     :param domain: domain of the `Data` object(s). If not specified, the
275     domain of the given `Data` objects is used.
276     :type domain: `escript.Domain`
277     :param timeStep: current timestep or sequence number - first one must be 0
278     :type timeStep: `int`
279     :param deltaT: timestep or sequence increment, see example above
280     :type deltaT: `int`
281     :param dynamicMesh: by default the mesh is assumed to be static and thus
282 caltinay 2599 only saved once at timestep 0 to save disk space.
283     Setting this to 1 changes the behaviour and the mesh
284     is saved at each timestep.
285 jfenwick 2625 :type dynamicMesh: `int`
286     :keyword <name>: writes the assigned value to the file using <name> as
287 caltinay 2188 identifier
288 jfenwick 2625 :type <name>: `Data` object
289 caltinay 2938 :note: The ESD concept is experimental and the file format likely to
290     change so use this function with caution.
291     :note: The data objects have to be defined on the same domain (but not
292     necessarily on the same `FunctionSpace`).
293 jfenwick 2625 :note: When saving a time series the first timestep must be 0 and it is
294 caltinay 2526 assumed that data from all timesteps share the domain. The dataset
295     file is updated in each iteration.
296 caltinay 2188 """
297     new_data = {}
298     for n,d in data.items():
299     if not d.isEmpty():
300     fs = d.getFunctionSpace()
301     domain2 = fs.getDomain()
302     if fs == escript.Solution(domain2):
303     new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
304     elif fs == escript.ReducedSolution(domain2):
305     new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
306     else:
307     new_data[n]=d
308     if domain==None: domain=domain2
309     if domain==None:
310     raise ValueError, "saveESD: no domain detected."
311    
312     if domain.onMasterProcessor() and not os.path.isdir(dataDir):
313     os.mkdir(dataDir)
314    
315 caltinay 2941 meshFile = datasetName+"_mesh"
316 caltinay 2599 fileNumber = timeStep / deltaT
317 caltinay 2526
318 caltinay 2599 if dynamicMesh == 0:
319     # later timesteps reuse mesh from t=0
320     if timeStep == 0:
321 caltinay 2941 domain.dump(os.path.join(dataDir, meshFile + ".nc"))
322 caltinay 2599 else:
323     meshFile += ".%04d"
324 caltinay 2941 domain.dump(os.path.join(dataDir, (meshFile + ".nc") % fileNumber))
325 caltinay 2526
326 caltinay 2188 outputString = ""
327 caltinay 2526
328 caltinay 2188 if domain.onMasterProcessor():
329     outputString += "#escript datafile V1.0\n"
330     # number of timesteps (currently only 1 is supported)
331 caltinay 2526 outputString += "T=%d\n" % (fileNumber+1)
332     # timestep increment
333     outputString += "DT=1\n"
334 caltinay 2188 # name of the mesh file
335     outputString += "M=%s\n" % meshFile
336     # number of blocks (MPI size)
337     outputString += "N=%d\n" % domain.getMPISize()
338    
339     # now add the variables
340     for varName, d in new_data.items():
341 caltinay 2941 varFile = datasetName+"_"+varName+".%04d"
342     d.dump(os.path.join(dataDir, (varFile + ".nc") % fileNumber))
343 caltinay 2188 if domain.onMasterProcessor():
344     outputString += "V=%s:%s\n" % (varFile, varName)
345    
346     if domain.onMasterProcessor():
347 caltinay 2938 esdfile = open(os.path.join(dataDir, datasetName+".esd"), "w")
348 caltinay 2188 esdfile.write(outputString)
349     esdfile.close()
350    
351 gross 290 def kronecker(d=3):
352     """
353 jfenwick 2625 Returns the kronecker delta-symbol.
354 jgs 123
355 jfenwick 2625 :param d: dimension or an object that has the ``getDim`` method defining the
356 caltinay 2169 dimension
357 jfenwick 2625 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
358     :return: the object u of rank 2 with *u[i,j]=1* for *i=j* and *u[i,j]=0*
359 caltinay 2169 otherwise
360 jfenwick 2625 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 2
361 gross 290 """
362     return identityTensor(d)
363    
364     def identity(shape=()):
365     """
366 jfenwick 2625 Returns the ``shape`` x ``shape`` identity tensor.
367 gross 290
368 jfenwick 2625 :param shape: input shape for the identity tensor
369     :type shape: ``tuple`` of ``int``
370     :return: array whose shape is shape x shape where *u[i,k]=1* for *i=k* and
371     *u[i,k]=0* otherwise for len(shape)=1. If len(shape)=2:
372     *u[i,j,k,l]=1* for *i=k and j=l* and *u[i,j,k,l]=0* otherwise.
373     :rtype: ``numpy.ndarray`` of rank 1, rank 2 or rank 4
374     :raise ValueError: if len(shape)>2
375 gross 290 """
376     if len(shape)>0:
377 jfenwick 2455 out=numpy.zeros(shape+shape,numpy.float64)
378 gross 290 if len(shape)==1:
379     for i0 in range(shape[0]):
380     out[i0,i0]=1.
381     elif len(shape)==2:
382     for i0 in range(shape[0]):
383     for i1 in range(shape[1]):
384     out[i0,i1,i0,i1]=1.
385     else:
386     raise ValueError,"identity: length of shape is restricted to 2."
387     else:
388     out=1.
389     return out
390    
391 gross 2912 def zeros(shape=()):
392     """
393     Returns the ``shape`` zero tensor.
394    
395     :param shape: input shape for the identity tensor
396     :type shape: ``tuple`` of ``int``
397     :return: array of shape filled with zeros
398     :rtype: ``numpy.ndarray``
399     """
400     if len(shape)>0:
401     out=numpy.zeros(shape,numpy.float64)
402     else:
403     out=0.
404     return out
405    
406 gross 290 def identityTensor(d=3):
407     """
408 jfenwick 2625 Returns the ``d`` x ``d`` identity matrix.
409 gross 290
410 jfenwick 2625 :param d: dimension or an object that has the ``getDim`` method defining the
411 caltinay 2169 dimension
412 jfenwick 2625 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
413     :return: the object u of rank 2 with *u[i,j]=1* for *i=j* and *u[i,j]=0*
414 caltinay 2169 otherwise
415 jfenwick 2625 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 2
416 gross 290 """
417 gross 442 if isinstance(d,escript.FunctionSpace):
418     return escript.Data(identity((d.getDim(),)),d)
419     elif isinstance(d,escript.Domain):
420     return identity((d.getDim(),))
421     else:
422     return identity((d,))
423 gross 290
424     def identityTensor4(d=3):
425     """
426 jfenwick 2625 Returns the ``d`` x ``d`` x ``d`` x ``d`` identity tensor.
427 gross 290
428 jfenwick 2625 :param d: dimension or an object that has the ``getDim`` method defining the
429 caltinay 2169 dimension
430 jfenwick 2625 :type d: ``int`` or any object with a ``getDim`` method
431     :return: the object u of rank 4 with *u[i,j,k,l]=1* for *i=k and j=l* and
432     *u[i,j,k,l]=0* otherwise
433     :rtype: ``numpy.ndarray`` or `escript.Data` of rank 4
434 gross 290 """
435 gross 442 if isinstance(d,escript.FunctionSpace):
436     return escript.Data(identity((d.getDim(),d.getDim())),d)
437     elif isinstance(d,escript.Domain):
438     return identity((d.getDim(),d.getDim()))
439     else:
440     return identity((d,d))
441 gross 290
442     def unitVector(i=0,d=3):
443     """
444 caltinay 2169 Returns a unit vector u of dimension d whose non-zero element is at index i.
445 gross 290
446 jfenwick 2625 :param i: index for non-zero element
447     :type i: ``int``
448     :param d: dimension or an object that has the ``getDim`` method defining the
449 caltinay 2169 dimension
450 jfenwick 2625 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
451     :return: the object u of rank 1 with *u[j]=1* for *j=index* and *u[j]=0*
452 caltinay 2169 otherwise
453 jfenwick 2625 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 1
454 gross 290 """
455     return kronecker(d)[i]
456    
457     #=========================================================================
458 caltinay 3493 # global reduction operations
459 gross 290 #=========================================================================
460     def Lsup(arg):
461 jgs 82 """
462 jfenwick 2625 Returns the Lsup-norm of argument ``arg``. This is the maximum absolute value
463     over all data points. This function is equivalent to ``sup(abs(arg))``.
464 jgs 149
465 jfenwick 2625 :param arg: argument
466     :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
467     :return: maximum value of the absolute value of ``arg`` over all components
468 caltinay 2169 and all data points
469 jfenwick 2625 :rtype: ``float``
470     :raise TypeError: if type of ``arg`` cannot be processed
471 jgs 123 """
472 jfenwick 2455 if isinstance(arg,numpy.ndarray):
473 gross 290 return sup(abs(arg))
474     elif isinstance(arg,escript.Data):
475     return arg._Lsup()
476     elif isinstance(arg,float):
477     return abs(arg)
478     elif isinstance(arg,int):
479     return abs(float(arg))
480 jgs 108 else:
481 caltinay 3496 raise TypeError,"Lsup: Unknown argument type."
482 jgs 82
483 gross 290 def sup(arg):
484 jgs 82 """
485 caltinay 2169 Returns the maximum value over all data points.
486 jgs 149
487 jfenwick 2625 :param arg: argument
488     :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
489     :return: maximum value of ``arg`` over all components and all data points
490     :rtype: ``float``
491     :raise TypeError: if type of ``arg`` cannot be processed
492 jgs 123 """
493 jfenwick 2455 if isinstance(arg,numpy.ndarray):
494 gross 290 return arg.max()
495     elif isinstance(arg,escript.Data):
496     return arg._sup()
497     elif isinstance(arg,float):
498     return arg
499     elif isinstance(arg,int):
500     return float(arg)
501 jgs 123 else:
502 caltinay 3496 raise TypeError,"sup: Unknown argument type."
503 jgs 82
504 gross 290 def inf(arg):
505 jgs 82 """
506 caltinay 2169 Returns the minimum value over all data points.
507 jgs 149
508 jfenwick 2625 :param arg: argument
509     :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
510     :return: minimum value of ``arg`` over all components and all data points
511     :rtype: ``float``
512     :raise TypeError: if type of ``arg`` cannot be processed
513 jgs 123 """
514 jfenwick 2455 if isinstance(arg,numpy.ndarray):
515 gross 290 return arg.min()
516     elif isinstance(arg,escript.Data):
517     return arg._inf()
518     elif isinstance(arg,float):
519     return arg
520     elif isinstance(arg,int):
521     return float(arg)
522 jgs 108 else:
523 gross 290 raise TypeError,"inf: Unknown argument type."
524 jgs 82
525 gross 290
526     #=========================================================================
527     # some little helpers
528     #=========================================================================
529 gross 912 def getRank(arg):
530 jgs 124 """
531 caltinay 2169 Identifies the rank of the argument.
532 gross 912
533 jfenwick 2625 :param arg: an object whose rank is to be returned
534 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``,
535     ``Symbol``
536 jfenwick 2625 :return: the rank of the argument
537     :rtype: ``int``
538     :raise TypeError: if type of ``arg`` cannot be processed
539 gross 912 """
540    
541 jfenwick 2455 if isinstance(arg,numpy.ndarray):
542     return arg.ndim
543 gross 912 elif isinstance(arg,escript.Data):
544     return arg.getRank()
545     elif isinstance(arg,float):
546     return 0
547     elif isinstance(arg,int):
548     return 0
549 caltinay 3507 elif isinstance(arg,Symbol):
550     return arg.getRank()
551 gross 912 else:
552 caltinay 2169 raise TypeError,"getRank: Unknown argument type."
553    
554 gross 912 def getShape(arg):
555     """
556 caltinay 2169 Identifies the shape of the argument.
557 jgs 149
558 jfenwick 2625 :param arg: an object whose shape is to be returned
559 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``,
560     ``Symbol``
561 jfenwick 2625 :return: the shape of the argument
562     :rtype: ``tuple`` of ``int``
563     :raise TypeError: if type of ``arg`` cannot be processed
564 jgs 124 """
565 gross 290
566 jfenwick 2455 if isinstance(arg,numpy.ndarray):
567 gross 290 return arg.shape
568 gross 3283 elif isinstance(arg,list):
569     return numpy.array(arg).shape
570 gross 290 elif isinstance(arg,escript.Data):
571     return arg.getShape()
572     elif isinstance(arg,float):
573     return ()
574     elif isinstance(arg,int):
575     return ()
576 caltinay 3507 elif isinstance(arg,Symbol):
577     return arg.getShape()
578 jgs 124 else:
579 caltinay 2169 raise TypeError,"getShape: Cannot identify shape"
580 jgs 124
581 gross 290 def pokeDim(arg):
582 jgs 82 """
583 caltinay 2169 Identifies the spatial dimension of the argument.
584 jgs 149
585 jfenwick 2625 :param arg: an object whose spatial dimension is to be returned
586     :type arg: any
587     :return: the spatial dimension of the argument, if available, or ``None``
588     :rtype: ``int`` or ``None``
589 jgs 123 """
590 gross 290
591     if isinstance(arg,escript.Data):
592     return arg.getFunctionSpace().getDim()
593 jgs 123 else:
594 gross 290 return None
595 jgs 82
596 caltinay 2169 def commonShape(arg0, arg1):
597 jgs 82 """
598 jfenwick 2625 Returns a shape to which ``arg0`` can be extended from the right and ``arg1``
599 caltinay 2169 can be extended from the left.
600 jgs 149
601 jfenwick 2625 :param arg0: an object with a shape (see `getShape`)
602     :param arg1: an object with a shape (see `getShape`)
603     :return: the shape of ``arg0`` or ``arg1`` such that the left part equals the
604     shape of ``arg0`` and the right end equals the shape of ``arg1``
605     :rtype: ``tuple`` of ``int``
606     :raise ValueError: if no shape can be found
607 jgs 123 """
608 gross 912 sh0=getShape(arg0)
609     sh1=getShape(arg1)
610 gross 290 if len(sh0)<len(sh1):
611     if not sh0==sh1[:len(sh0)]:
612     raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
613     return sh1
614     elif len(sh0)>len(sh1):
615     if not sh1==sh0[:len(sh1)]:
616     raise ValueError,"argument 1 cannot be extended to the shape of argument 0"
617     return sh0
618 jgs 82 else:
619 gross 290 if not sh0==sh1:
620     raise ValueError,"argument 1 and argument 0 have not the same shape."
621     return sh0
622 jgs 82
623 gross 290 def commonDim(*args):
624 jgs 82 """
625 caltinay 2169 Identifies, if possible, the spatial dimension across a set of objects
626     which may or may not have a spatial dimension.
627 jgs 149
628 jfenwick 2625 :param args: given objects
629     :return: the spatial dimension of the objects with identifiable dimension
630     (see `pokeDim`). If none of the objects has a spatial dimension
631     ``None`` is returned.
632     :rtype: ``int`` or ``None``
633     :raise ValueError: if the objects with identifiable dimension don't have
634 caltinay 2169 the same spatial dimension.
635 jgs 123 """
636 gross 290 out=None
637     for a in args:
638     d=pokeDim(a)
639     if not out==None:
640     if not (d==None or out==d):
641     raise ValueError,"dimension of arguments don't match"
642     else:
643     out=d
644     return out
645 jgs 82
646 gross 290 def testForZero(arg):
647 jgs 150 """
648 caltinay 2169 Tests if the argument is identical to zero.
649 jgs 123
650 jfenwick 2625 :param arg: the object to test for zero
651     :type arg: typically ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
652     :return: True if the argument is identical to zero, False otherwise
653     :rtype: ``bool``
654 jgs 150 """
655 jfenwick 2455 if isinstance(arg,numpy.ndarray):
656 gross 290 return not Lsup(arg)>0.
657 gross 396 elif isinstance(arg,escript.Data):
658 gross 290 return False
659 gross 396 elif isinstance(arg,float):
660     return not Lsup(arg)>0.
661     elif isinstance(arg,int):
662     return not Lsup(arg)>0.
663     else:
664     return False
665 jgs 150
666 gross 290 def matchType(arg0=0.,arg1=0.):
667 jgs 150 """
668 jfenwick 2625 Converts ``arg0`` and ``arg1`` both to the same type ``numpy.ndarray`` or
669 caltinay 3493 `escript.Data`
670 jgs 150
671 jfenwick 2625 :param arg0: first argument
672 caltinay 3507 :type arg0: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, ``Symbol``
673 jfenwick 2625 :param arg1: second argument
674 caltinay 3507 :type arg1: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, ``Symbol``
675     :return: a tuple representing ``arg0`` and ``arg1`` with the same type or
676     with at least one of them being a `Symbol`
677 caltinay 3493 :rtype: ``tuple`` of two ``numpy.ndarray`` or two `escript.Data`
678 jfenwick 2625 :raise TypeError: if type of ``arg0`` or ``arg1`` cannot be processed
679 jgs 150 """
680 jfenwick 2455 if isinstance(arg0,numpy.ndarray):
681     if isinstance(arg1,numpy.ndarray):
682 gross 290 pass
683     elif isinstance(arg1,escript.Data):
684     arg0=escript.Data(arg0,arg1.getFunctionSpace())
685     elif isinstance(arg1,float):
686 jfenwick 2455 arg1=numpy.array(arg1,dtype=numpy.float64)
687 gross 290 elif isinstance(arg1,int):
688 jfenwick 2455 arg1=numpy.array(float(arg1),dtype=numpy.float64)
689 caltinay 3507 elif isinstance(arg1,Symbol):
690     pass
691 gross 290 else:
692 caltinay 2169 raise TypeError,"function: Unknown type of second argument."
693 gross 290 elif isinstance(arg0,escript.Data):
694 jfenwick 2455 if isinstance(arg1,numpy.ndarray):
695 gross 290 arg1=escript.Data(arg1,arg0.getFunctionSpace())
696     elif isinstance(arg1,escript.Data):
697     pass
698     elif isinstance(arg1,float):
699     arg1=escript.Data(arg1,(),arg0.getFunctionSpace())
700     elif isinstance(arg1,int):
701     arg1=escript.Data(float(arg1),(),arg0.getFunctionSpace())
702 caltinay 3507 elif isinstance(arg1,Symbol):
703     pass
704 gross 290 else:
705 caltinay 2169 raise TypeError,"function: Unknown type of second argument."
706 caltinay 3507 elif isinstance(arg0,Symbol):
707     if isinstance(arg1,numpy.ndarray):
708     pass
709     elif isinstance(arg1,escript.Data):
710     pass
711     elif isinstance(arg1,float):
712     pass
713     elif isinstance(arg1,int):
714     pass
715     elif isinstance(arg1,Symbol):
716     pass
717     else:
718     raise TypeError,"function: Unknown type of second argument."
719 gross 290 elif isinstance(arg0,float):
720 jfenwick 2455 if isinstance(arg1,numpy.ndarray):
721     arg0=numpy.array(arg0,dtype=numpy.float64)
722 gross 290 elif isinstance(arg1,escript.Data):
723     arg0=escript.Data(arg0,arg1.getFunctionSpace())
724     elif isinstance(arg1,float):
725 jfenwick 2455 arg0=numpy.array(arg0,dtype=numpy.float64)
726     arg1=numpy.array(arg1,dtype=numpy.float64)
727 gross 290 elif isinstance(arg1,int):
728 jfenwick 2455 arg0=numpy.array(arg0,dtype=numpy.float64)
729     arg1=numpy.array(float(arg1),dtype=numpy.float64)
730 caltinay 3507 elif isinstance(arg1,Symbol):
731     pass
732 gross 290 else:
733 caltinay 2169 raise TypeError,"function: Unknown type of second argument."
734 gross 290 elif isinstance(arg0,int):
735 jfenwick 2455 if isinstance(arg1,numpy.ndarray):
736     arg0=numpy.array(float(arg0),dtype=numpy.float64)
737 gross 290 elif isinstance(arg1,escript.Data):
738     arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
739     elif isinstance(arg1,float):
740 jfenwick 2455 arg0=numpy.array(float(arg0),dtype=numpy.float64)
741     arg1=numpy.array(arg1,dtype=numpy.float64)
742 gross 290 elif isinstance(arg1,int):
743 jfenwick 2455 arg0=numpy.array(float(arg0),dtype=numpy.float64)
744     arg1=numpy.array(float(arg1),dtype=numpy.float64)
745 caltinay 3507 elif isinstance(arg1,Symbol):
746     pass
747 gross 290 else:
748 caltinay 2169 raise TypeError,"function: Unknown type of second argument."
749 jgs 150 else:
750 caltinay 2169 raise TypeError,"function: Unknown type of first argument."
751 jgs 150
752 gross 290 return arg0,arg1
753    
754     def matchShape(arg0,arg1):
755 jgs 150 """
756 jfenwick 2625 Returns a representation of ``arg0`` and ``arg1`` which have the same shape.
757 jgs 150
758 jfenwick 2625 :param arg0: first argument
759 caltinay 3507 :type arg0: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, `Symbol`
760 jfenwick 2625 :param arg1: second argument
761 caltinay 3507 :type arg1: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, `Symbol`
762 jfenwick 2625 :return: ``arg0`` and ``arg1`` where copies are returned when the shape has
763 caltinay 2169 to be changed
764 jfenwick 2625 :rtype: ``tuple``
765 jgs 150 """
766 gross 290 sh=commonShape(arg0,arg1)
767 gross 912 sh0=getShape(arg0)
768     sh1=getShape(arg1)
769 gross 290 if len(sh0)<len(sh):
770 jfenwick 2455 return outer(arg0,numpy.ones(sh[len(sh0):],numpy.float64)),arg1
771 gross 290 elif len(sh1)<len(sh):
772 jfenwick 2455 return arg0,outer(arg1,numpy.ones(sh[len(sh1):],numpy.float64))
773 caltinay 2169 else:
774 gross 290 return arg0,arg1
775 caltinay 2169
776 gross 290 def log10(arg):
777     """
778 jfenwick 2625 Returns base-10 logarithm of argument ``arg``.
779 gross 290
780 jfenwick 2625 :param arg: argument
781 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
782     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
783 jfenwick 2625 on the type of ``arg``
784     :raise TypeError: if the type of the argument is not expected
785 gross 290 """
786 jfenwick 2455 if isinstance(arg,numpy.ndarray):
787     return numpy.log10(arg)
788 gross 290 elif isinstance(arg,escript.Data):
789     return arg._log10()
790     elif isinstance(arg,float):
791     return math.log10(arg)
792     elif isinstance(arg,int):
793     return math.log10(float(arg))
794 caltinay 3507 elif isinstance(arg,Symbol):
795     return arg.applyfunc(symfn.log10)
796 gross 290 else:
797     raise TypeError,"log10: Unknown argument type."
798    
799     def wherePositive(arg):
800     """
801 jfenwick 2625 Returns mask of positive values of argument ``arg``.
802 gross 290
803 jfenwick 2625 :param arg: argument
804 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
805     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
806 jfenwick 2625 on the type of ``arg``
807     :raise TypeError: if the type of the argument is not expected
808 gross 290 """
809 jfenwick 2455 if isinstance(arg,numpy.ndarray):
810     out=numpy.greater(arg,numpy.zeros(arg.shape,numpy.float64))*1.
811     if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
812 gross 396 return out
813 gross 290 elif isinstance(arg,escript.Data):
814     return arg._wherePositive()
815     elif isinstance(arg,float):
816     if arg>0:
817     return 1.
818     else:
819     return 0.
820     elif isinstance(arg,int):
821     if arg>0:
822     return 1.
823     else:
824     return 0.
825 caltinay 3507 elif isinstance(arg,Symbol):
826     return arg.applyfunc(symfn.wherePositive)
827 gross 290 else:
828     raise TypeError,"wherePositive: Unknown argument type."
829    
830     def whereNegative(arg):
831     """
832 jfenwick 2625 Returns mask of negative values of argument ``arg``.
833 gross 290
834 jfenwick 2625 :param arg: argument
835 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
836     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
837 jfenwick 2625 on the type of ``arg``
838     :raise TypeError: if the type of the argument is not expected
839 gross 290 """
840 jfenwick 2455 if isinstance(arg,numpy.ndarray):
841     out=numpy.less(arg,numpy.zeros(arg.shape,numpy.float64))*1.
842     if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
843 gross 396 return out
844 gross 290 elif isinstance(arg,escript.Data):
845     return arg._whereNegative()
846     elif isinstance(arg,float):
847     if arg<0:
848     return 1.
849     else:
850     return 0.
851     elif isinstance(arg,int):
852     if arg<0:
853     return 1.
854     else:
855     return 0.
856 caltinay 3507 elif isinstance(arg,Symbol):
857     return arg.applyfunc(symfn.whereNegative)
858 gross 290 else:
859     raise TypeError,"whereNegative: Unknown argument type."
860    
861     def whereNonNegative(arg):
862     """
863 jfenwick 2625 Returns mask of non-negative values of argument ``arg``.
864 gross 290
865 jfenwick 2625 :param arg: argument
866 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
867     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
868 jfenwick 2625 on the type of ``arg``
869     :raise TypeError: if the type of the argument is not expected
870 gross 290 """
871 jfenwick 2455 if isinstance(arg,numpy.ndarray):
872     out=numpy.greater_equal(arg,numpy.zeros(arg.shape,numpy.float64))*1.
873     if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
874 gross 396 return out
875 gross 290 elif isinstance(arg,escript.Data):
876     return arg._whereNonNegative()
877     elif isinstance(arg,float):
878     if arg<0:
879     return 0.
880     else:
881     return 1.
882     elif isinstance(arg,int):
883     if arg<0:
884     return 0.
885     else:
886     return 1.
887 caltinay 3507 elif isinstance(arg,Symbol):
888     return arg.applyfunc(symfn.whereNonNegative)
889 gross 290 else:
890     raise TypeError,"whereNonNegative: Unknown argument type."
891    
892     def whereNonPositive(arg):
893     """
894 jfenwick 2625 Returns mask of non-positive values of argument ``arg``.
895 gross 290
896 jfenwick 2625 :param arg: argument
897 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
898     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
899 jfenwick 2625 on the type of ``arg``
900     :raise TypeError: if the type of the argument is not expected
901 gross 290 """
902 jfenwick 2455 if isinstance(arg,numpy.ndarray):
903     out=numpy.less_equal(arg,numpy.zeros(arg.shape,numpy.float64))*1.
904     if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
905 gross 396 return out
906 gross 290 elif isinstance(arg,escript.Data):
907     return arg._whereNonPositive()
908     elif isinstance(arg,float):
909     if arg>0:
910     return 0.
911     else:
912     return 1.
913     elif isinstance(arg,int):
914     if arg>0:
915     return 0.
916     else:
917     return 1.
918 caltinay 3507 elif isinstance(arg,Symbol):
919     return arg.applyfunc(symfn.whereNonPositive)
920 gross 290 else:
921     raise TypeError,"whereNonPositive: Unknown argument type."
922    
923 gross 3358 def whereZero(arg,tol=None,rtol=math.sqrt(EPSILON)):
924 gross 290 """
925 jfenwick 2625 Returns mask of zero entries of argument ``arg``.
926 gross 290
927 jfenwick 2625 :param arg: argument
928 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
929 jfenwick 2625 :param tol: absolute tolerance. Values with absolute value less than tol are accepted
930     as zero. If ``tol`` is not present ``rtol``*```Lsup` (arg)`` is used.
931     :type tol: ``float``
932     :param rtol: relative tolerance used to define the absolute tolerance if ``tol`` is not present.
933     :type rtol: non-negative ``float``
934 caltinay 3507 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
935 jfenwick 2625 on the type of ``arg``
936     :raise ValueError: if ``rtol`` is non-negative.
937     :raise TypeError: if the type of the argument is not expected
938 gross 290 """
939 caltinay 3507 if tol == None and not isinstance(arg, Symbol):
940 caltinay 3493 if rtol<0: raise ValueError,"rtol must be non-negative."
941     tol = Lsup(arg)*rtol
942 jfenwick 2455 if isinstance(arg,numpy.ndarray):
943     out=numpy.less_equal(abs(arg)-tol,numpy.zeros(arg.shape,numpy.float64))*1.
944     if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
945 gross 396 return out
946 gross 290 elif isinstance(arg,escript.Data):
947 gross 698 return arg._whereZero(tol)
948 gross 290 elif isinstance(arg,float):
949     if abs(arg)<=tol:
950     return 1.
951     else:
952     return 0.
953     elif isinstance(arg,int):
954     if abs(float(arg))<=tol:
955     return 1.
956     else:
957     return 0.
958 caltinay 3507 elif isinstance(arg,Symbol):
959     return arg.applyfunc(symfn.whereZero)
960 gross 290 else:
961     raise TypeError,"whereZero: Unknown argument type."
962    
963     def whereNonZero(arg,tol=0.):
964     """
965 jfenwick 2625 Returns mask of values different from zero of argument ``arg``.
966 gross 290
967 jfenwick 2625 :param arg: argument
968 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
969 jfenwick 2625 :param tol: absolute tolerance. Values with absolute value less than tol are accepted
970     as zero. If ``tol`` is not present ``rtol``*```Lsup` (arg)`` is used.
971     :type tol: ``float``
972 caltinay 3507 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
973 jfenwick 2625 on the type of ``arg``
974     :raise ValueError: if ``rtol`` is non-negative.
975     :raise TypeError: if the type of the argument is not expected
976 gross 290 """
977 gross 2304 if tol == None:
978 caltinay 3493 if rtol<=0: raise ValueError,"rtol must be non-negative."
979     tol = Lsup(arg)*rtol
980 jfenwick 2455 if isinstance(arg,numpy.ndarray):
981     out=numpy.greater(abs(arg)-tol,numpy.zeros(arg.shape,numpy.float64))*1.
982     if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
983 gross 396 return out
984 gross 290 elif isinstance(arg,escript.Data):
985 gross 698 return arg._whereNonZero(tol)
986 gross 290 elif isinstance(arg,float):
987     if abs(arg)>tol:
988     return 1.
989     else:
990     return 0.
991     elif isinstance(arg,int):
992     if abs(float(arg))>tol:
993     return 1.
994     else:
995     return 0.
996 caltinay 3507 elif isinstance(arg,Symbol):
997     return arg.applyfunc(symfn.whereNonZero)
998 gross 290 else:
999     raise TypeError,"whereNonZero: Unknown argument type."
1000    
1001 ksteube 876 def erf(arg):
1002     """
1003 jfenwick 2625 Returns the error function *erf* of argument ``arg``.
1004 ksteube 876
1005 jfenwick 2625 :param arg: argument
1006 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1007     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1008 jfenwick 2625 on the type of ``arg``
1009     :raise TypeError: if the type of the argument is not expected
1010 ksteube 876 """
1011     if isinstance(arg,escript.Data):
1012     return arg._erf()
1013 caltinay 3507 elif isinstance(arg,Symbol):
1014     return arg.applyfunc(symfn.erf)
1015 ksteube 876 else:
1016     raise TypeError,"erf: Unknown argument type."
1017    
1018 gross 290 def sin(arg):
1019     """
1020 jfenwick 2625 Returns sine of argument ``arg``.
1021 gross 290
1022 jfenwick 2625 :param arg: argument
1023 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1024     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1025 jfenwick 2625 on the type of ``arg``
1026     :raise TypeError: if the type of the argument is not expected
1027 gross 290 """
1028 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1029     return numpy.sin(arg)
1030 gross 290 elif isinstance(arg,escript.Data):
1031     return arg._sin()
1032     elif isinstance(arg,float):
1033     return math.sin(arg)
1034     elif isinstance(arg,int):
1035     return math.sin(arg)
1036 caltinay 3507 elif isinstance(arg,Symbol):
1037     return arg.applyfunc(symfn.sin)
1038 gross 290 else:
1039     raise TypeError,"sin: Unknown argument type."
1040    
1041     def cos(arg):
1042     """
1043 jfenwick 2625 Returns cosine of argument ``arg``.
1044 gross 290
1045 jfenwick 2625 :param arg: argument
1046 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1047     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1048 jfenwick 2625 on the type of ``arg``
1049     :raise TypeError: if the type of the argument is not expected
1050 gross 290 """
1051 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1052     return numpy.cos(arg)
1053 gross 290 elif isinstance(arg,escript.Data):
1054     return arg._cos()
1055     elif isinstance(arg,float):
1056     return math.cos(arg)
1057     elif isinstance(arg,int):
1058     return math.cos(arg)
1059 caltinay 3507 elif isinstance(arg,Symbol):
1060     return arg.applyfunc(symfn.cos)
1061 gross 290 else:
1062     raise TypeError,"cos: Unknown argument type."
1063    
1064     def tan(arg):
1065     """
1066 jfenwick 2625 Returns tangent of argument ``arg``.
1067 gross 290
1068 jfenwick 2625 :param arg: argument
1069 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1070     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1071 jfenwick 2625 on the type of ``arg``
1072     :raise TypeError: if the type of the argument is not expected
1073 gross 290 """
1074 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1075     return numpy.tan(arg)
1076 gross 290 elif isinstance(arg,escript.Data):
1077     return arg._tan()
1078     elif isinstance(arg,float):
1079     return math.tan(arg)
1080     elif isinstance(arg,int):
1081     return math.tan(arg)
1082 caltinay 3507 elif isinstance(arg,Symbol):
1083     return arg.applyfunc(symfn.tan)
1084 gross 290 else:
1085     raise TypeError,"tan: Unknown argument type."
1086    
1087     def asin(arg):
1088     """
1089 jfenwick 2625 Returns the inverse sine of argument ``arg``.
1090 gross 290
1091 jfenwick 2625 :param arg: argument
1092 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1093     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1094 jfenwick 2625 on the type of ``arg``
1095     :raise TypeError: if the type of the argument is not expected
1096 gross 290 """
1097 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1098     return numpy.arcsin(arg)
1099 gross 290 elif isinstance(arg,escript.Data):
1100     return arg._asin()
1101     elif isinstance(arg,float):
1102     return math.asin(arg)
1103     elif isinstance(arg,int):
1104     return math.asin(arg)
1105 caltinay 3507 elif isinstance(arg,Symbol):
1106     return arg.applyfunc(symfn.asin)
1107 gross 290 else:
1108     raise TypeError,"asin: Unknown argument type."
1109    
1110     def acos(arg):
1111     """
1112 jfenwick 2625 Returns the inverse cosine of argument ``arg``.
1113 gross 290
1114 jfenwick 2625 :param arg: argument
1115 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1116     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1117 jfenwick 2625 on the type of ``arg``
1118     :raise TypeError: if the type of the argument is not expected
1119 gross 290 """
1120 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1121     return numpy.arccos(arg)
1122 gross 290 elif isinstance(arg,escript.Data):
1123     return arg._acos()
1124     elif isinstance(arg,float):
1125     return math.acos(arg)
1126     elif isinstance(arg,int):
1127     return math.acos(arg)
1128 caltinay 3507 elif isinstance(arg,Symbol):
1129     return arg.applyfunc(symfn.acos)
1130 gross 290 else:
1131     raise TypeError,"acos: Unknown argument type."
1132    
1133     def atan(arg):
1134     """
1135 jfenwick 2625 Returns inverse tangent of argument ``arg``.
1136 gross 290
1137 jfenwick 2625 :param arg: argument
1138 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1139     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1140 jfenwick 2625 on the type of ``arg``
1141     :raise TypeError: if the type of the argument is not expected
1142 gross 290 """
1143 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1144     return numpy.arctan(arg)
1145 gross 290 elif isinstance(arg,escript.Data):
1146     return arg._atan()
1147     elif isinstance(arg,float):
1148     return math.atan(arg)
1149     elif isinstance(arg,int):
1150     return math.atan(arg)
1151 caltinay 3507 elif isinstance(arg,Symbol):
1152     return arg.applyfunc(symfn.atan)
1153 gross 290 else:
1154     raise TypeError,"atan: Unknown argument type."
1155    
1156 jgs 150 def sinh(arg):
1157 gross 290 """
1158 jfenwick 2625 Returns the hyperbolic sine of argument ``arg``.
1159 jgs 150
1160 jfenwick 2625 :param arg: argument
1161 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1162     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1163 jfenwick 2625 on the type of ``arg``
1164     :raise TypeError: if the type of the argument is not expected
1165 gross 290 """
1166 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1167     return numpy.sinh(arg)
1168 gross 290 elif isinstance(arg,escript.Data):
1169     return arg._sinh()
1170     elif isinstance(arg,float):
1171     return math.sinh(arg)
1172     elif isinstance(arg,int):
1173     return math.sinh(arg)
1174 caltinay 3507 elif isinstance(arg,Symbol):
1175     return arg.applyfunc(symfn.sinh)
1176 gross 290 else:
1177     raise TypeError,"sinh: Unknown argument type."
1178 jgs 150
1179     def cosh(arg):
1180 gross 290 """
1181 jfenwick 2625 Returns the hyperbolic cosine of argument ``arg``.
1182 jgs 150
1183 jfenwick 2625 :param arg: argument
1184 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1185     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1186 jfenwick 2625 on the type of ``arg``
1187     :raise TypeError: if the type of the argument is not expected
1188 gross 290 """
1189 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1190     return numpy.cosh(arg)
1191 gross 290 elif isinstance(arg,escript.Data):
1192     return arg._cosh()
1193     elif isinstance(arg,float):
1194     return math.cosh(arg)
1195     elif isinstance(arg,int):
1196     return math.cosh(arg)
1197 caltinay 3507 elif isinstance(arg,Symbol):
1198     return arg.applyfunc(symfn.cosh)
1199 gross 290 else:
1200     raise TypeError,"cosh: Unknown argument type."
1201 jgs 150
1202     def tanh(arg):
1203 gross 290 """
1204 jfenwick 2625 Returns the hyperbolic tangent of argument ``arg``.
1205 jgs 150
1206 jfenwick 2625 :param arg: argument
1207 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1208     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1209 jfenwick 2625 on the type of ``arg``
1210     :raise TypeError: if the type of the argument is not expected
1211 gross 290 """
1212 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1213     return numpy.tanh(arg)
1214 gross 290 elif isinstance(arg,escript.Data):
1215     return arg._tanh()
1216     elif isinstance(arg,float):
1217     return math.tanh(arg)
1218     elif isinstance(arg,int):
1219     return math.tanh(arg)
1220 caltinay 3507 elif isinstance(arg,Symbol):
1221     return arg.applyfunc(symfn.tanh)
1222 gross 290 else:
1223     raise TypeError,"tanh: Unknown argument type."
1224 jgs 150
1225     def asinh(arg):
1226 gross 290 """
1227 jfenwick 2625 Returns the inverse hyperbolic sine of argument ``arg``.
1228 jgs 150
1229 jfenwick 2625 :param arg: argument
1230 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1231     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1232 jfenwick 2625 on the type of ``arg``
1233     :raise TypeError: if the type of the argument is not expected
1234 gross 290 """
1235 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1236     return numpy.arcsinh(arg)
1237 gross 290 elif isinstance(arg,escript.Data):
1238     return arg._asinh()
1239     elif isinstance(arg,float):
1240 jfenwick 2455 return numpy.arcsinh(arg)
1241 gross 290 elif isinstance(arg,int):
1242 jfenwick 2455 return numpy.arcsinh(float(arg))
1243 caltinay 3507 elif isinstance(arg,Symbol):
1244     return arg.applyfunc(symfn.asinh)
1245 gross 290 else:
1246     raise TypeError,"asinh: Unknown argument type."
1247 jgs 150
1248     def acosh(arg):
1249 gross 290 """
1250 jfenwick 2625 Returns the inverse hyperbolic cosine of argument ``arg``.
1251 jgs 150
1252 jfenwick 2625 :param arg: argument
1253 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1254     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1255 jfenwick 2625 on the type of ``arg``
1256     :raise TypeError: if the type of the argument is not expected
1257 gross 290 """
1258 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1259     return numpy.arccosh(arg)
1260 gross 290 elif isinstance(arg,escript.Data):
1261     return arg._acosh()
1262     elif isinstance(arg,float):
1263 jfenwick 2455 return numpy.arccosh(arg)
1264 gross 290 elif isinstance(arg,int):
1265 jfenwick 2455 return numpy.arccosh(float(arg))
1266 caltinay 3507 elif isinstance(arg,Symbol):
1267     return arg.applyfunc(symfn.acosh)
1268 gross 290 else:
1269     raise TypeError,"acosh: Unknown argument type."
1270 jgs 150
1271     def atanh(arg):
1272 gross 290 """
1273 jfenwick 2625 Returns the inverse hyperbolic tangent of argument ``arg``.
1274 jgs 150
1275 jfenwick 2625 :param arg: argument
1276 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1277     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1278 jfenwick 2625 on the type of ``arg``
1279     :raise TypeError: if the type of the argument is not expected
1280 gross 290 """
1281 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1282     return numpy.arctanh(arg)
1283 gross 290 elif isinstance(arg,escript.Data):
1284     return arg._atanh()
1285     elif isinstance(arg,float):
1286 jfenwick 2455 return numpy.arctanh(arg)
1287 gross 290 elif isinstance(arg,int):
1288 jfenwick 2455 return numpy.arctanh(float(arg))
1289 caltinay 3507 elif isinstance(arg,Symbol):
1290     return arg.applyfunc(symfn.atanh)
1291 gross 290 else:
1292     raise TypeError,"atanh: Unknown argument type."
1293 jgs 150
1294 gross 290 def exp(arg):
1295     """
1296 jfenwick 2625 Returns *e* to the power of argument ``arg``.
1297 gross 290
1298 jfenwick 2625 :param arg: argument
1299 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1300     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1301 caltinay 2169 on the type of arg
1302 jfenwick 2625 :raise TypeError: if the type of the argument is not expected
1303 gross 290 """
1304 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1305     return numpy.exp(arg)
1306 gross 290 elif isinstance(arg,escript.Data):
1307     return arg._exp()
1308     elif isinstance(arg,float):
1309     return math.exp(arg)
1310     elif isinstance(arg,int):
1311     return math.exp(arg)
1312 caltinay 3507 elif isinstance(arg,Symbol):
1313     return arg.applyfunc(symfn.exp)
1314 gross 290 else:
1315     raise TypeError,"exp: Unknown argument type."
1316    
1317     def sqrt(arg):
1318     """
1319 jfenwick 2625 Returns the square root of argument ``arg``.
1320 gross 290
1321 jfenwick 2625 :param arg: argument
1322 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1323     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1324 jfenwick 2625 depending on the type of ``arg``
1325     :raise TypeError: if the type of the argument is not expected
1326 gross 290 """
1327 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1328     return numpy.sqrt(arg)
1329 gross 290 elif isinstance(arg,escript.Data):
1330     return arg._sqrt()
1331     elif isinstance(arg,float):
1332     return math.sqrt(arg)
1333     elif isinstance(arg,int):
1334     return math.sqrt(arg)
1335 caltinay 3507 elif isinstance(arg,Symbol):
1336     return arg.applyfunc(symfn.sqrt)
1337 gross 290 else:
1338     raise TypeError,"sqrt: Unknown argument type."
1339    
1340     def log(arg):
1341     """
1342 jfenwick 2625 Returns the natural logarithm of argument ``arg``.
1343 gross 290
1344 jfenwick 2625 :param arg: argument
1345 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1346     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1347 jfenwick 2625 on the type of ``arg``
1348     :raise TypeError: if the type of the argument is not expected
1349 gross 290 """
1350 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1351     return numpy.log(arg)
1352 gross 290 elif isinstance(arg,escript.Data):
1353     return arg._log()
1354     elif isinstance(arg,float):
1355     return math.log(arg)
1356     elif isinstance(arg,int):
1357     return math.log(arg)
1358 caltinay 3507 elif isinstance(arg,Symbol):
1359     return arg.applyfunc(symfn.log)
1360 gross 290 else:
1361     raise TypeError,"log: Unknown argument type."
1362    
1363 jgs 123 def sign(arg):
1364 gross 290 """
1365 jfenwick 2625 Returns the sign of argument ``arg``.
1366 jgs 149
1367 jfenwick 2625 :param arg: argument
1368 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1369     :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1370 jfenwick 2625 on the type of ``arg``
1371     :raise TypeError: if the type of the argument is not expected
1372 gross 290 """
1373 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1374 gross 329 return wherePositive(arg)-whereNegative(arg)
1375 gross 290 elif isinstance(arg,escript.Data):
1376     return arg._sign()
1377     elif isinstance(arg,float):
1378     if arg>0:
1379     return 1.
1380     elif arg<0:
1381     return -1.
1382     else:
1383     return 0.
1384     elif isinstance(arg,int):
1385     if float(arg)>0:
1386     return 1.
1387     elif float(arg)<0:
1388     return -1.
1389     else:
1390     return 0.
1391 caltinay 3507 elif isinstance(arg,Symbol):
1392     return arg.applyfunc(symfn.sign)
1393 gross 290 else:
1394     raise TypeError,"sign: Unknown argument type."
1395 jgs 82
1396 gross 312 def minval(arg):
1397     """
1398 jfenwick 2625 Returns the minimum value over all components of ``arg`` at each data point.
1399 gross 290
1400 jfenwick 2625 :param arg: argument
1401 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1402     :rtype: ``float``, `escript.Data`, `Symbol` depending on the type of ``arg``
1403 jfenwick 2625 :raise TypeError: if the type of the argument is not expected
1404 gross 312 """
1405 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1406     if arg.ndim==0:
1407 gross 341 return float(arg)
1408     else:
1409     return arg.min()
1410 gross 312 elif isinstance(arg,escript.Data):
1411     return arg._minval()
1412     elif isinstance(arg,float):
1413     return arg
1414     elif isinstance(arg,int):
1415     return float(arg)
1416 caltinay 3507 elif isinstance(arg,Symbol):
1417 caltinay 3493 return symfn.minval(arg)
1418 gross 312 else:
1419     raise TypeError,"minval: Unknown argument type."
1420    
1421     def maxval(arg):
1422     """
1423 jfenwick 2625 Returns the maximum value over all components of ``arg`` at each data point.
1424 gross 312
1425 jfenwick 2625 :param arg: argument
1426 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1427     :rtype: ``float``, `escript.Data`, `Symbol` depending on the type of ``arg``
1428 jfenwick 2625 :raise TypeError: if the type of the argument is not expected
1429 gross 312 """
1430 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1431     if arg.ndim==0:
1432 gross 341 return float(arg)
1433     else:
1434     return arg.max()
1435 gross 312 elif isinstance(arg,escript.Data):
1436     return arg._maxval()
1437     elif isinstance(arg,float):
1438     return arg
1439     elif isinstance(arg,int):
1440     return float(arg)
1441 caltinay 3507 elif isinstance(arg,Symbol):
1442 caltinay 3493 return symfn.maxval(arg)
1443 gross 312 else:
1444     raise TypeError,"maxval: Unknown argument type."
1445    
1446     def length(arg):
1447     """
1448 jfenwick 2625 Returns the length (Euclidean norm) of argument ``arg`` at each data point.
1449 gross 312
1450 jfenwick 2625 :param arg: argument
1451 caltinay 3507 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1452     :rtype: ``float``, `escript.Data`, `Symbol` depending on the type of ``arg``
1453 gross 312 """
1454     return sqrt(inner(arg,arg))
1455    
1456 gross 429 def trace(arg,axis_offset=0):
1457     """
1458 jfenwick 2625 Returns the trace of ``arg`` which is the sum of ``arg[k,k]`` over k.
1459 gross 429
1460 jfenwick 2625 :param arg: argument
1461 caltinay 3507 :type arg: `escript.Data`, `Symbol`, ``numpy.ndarray``
1462 jfenwick 2625 :param axis_offset: ``axis_offset`` to components to sum over. ``axis_offset``
1463     must be non-negative and less than the rank of ``arg`` +1.
1464     The dimensions of component ``axis_offset`` and
1465 caltinay 2169 axis_offset+1 must be equal.
1466 jfenwick 2625 :type axis_offset: ``int``
1467     :return: trace of arg. The rank of the returned object is rank of ``arg``
1468 caltinay 2169 minus 2.
1469 caltinay 3507 :rtype: `escript.Data`, `Symbol` or ``numpy.ndarray`` depending on the type
1470     of ``arg``
1471 gross 429 """
1472 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1473 gross 429 sh=arg.shape
1474 caltinay 2169 if len(sh)<2:
1475 gross 785 raise ValueError,"rank of argument must be greater than 1"
1476 gross 429 if axis_offset<0 or axis_offset>len(sh)-2:
1477 jfenwick 2200 raise ValueError,"axis_offset must be between 0 and %d"%(len(sh)-2)
1478 gross 429 s1=1
1479     for i in range(axis_offset): s1*=sh[i]
1480     s2=1
1481     for i in range(axis_offset+2,len(sh)): s2*=sh[i]
1482     if not sh[axis_offset] == sh[axis_offset+1]:
1483 jfenwick 2200 raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
1484 jfenwick 2455 arg_reshaped=numpy.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
1485     out=numpy.zeros([s1,s2],numpy.float64)
1486 gross 429 for i1 in range(s1):
1487     for i2 in range(s2):
1488     for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
1489     out.resize(sh[:axis_offset]+sh[axis_offset+2:])
1490     return out
1491     elif isinstance(arg,escript.Data):
1492 caltinay 2169 if arg.getRank()<2:
1493 gross 785 raise ValueError,"rank of argument must be greater than 1"
1494     if axis_offset<0 or axis_offset>arg.getRank()-2:
1495 jfenwick 2200 raise ValueError,"axis_offset must be between 0 and %d"%(arg.getRank()-2)
1496 caltinay 2169 s=list(arg.getShape())
1497 gross 785 if not s[axis_offset] == s[axis_offset+1]:
1498 jfenwick 2200 raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
1499 gross 800 return arg._trace(axis_offset)
1500 caltinay 3507 elif isinstance(arg,Symbol):
1501     if arg.getRank()<2:
1502     raise ValueError,"rank of argument must be greater than 1"
1503     if axis_offset<0 or axis_offset>arg.getRank()-2:
1504     raise ValueError,"axis_offset must be between 0 and %d"%(arg.getRank()-2)
1505     s=list(arg.getShape())
1506     if not s[axis_offset] == s[axis_offset+1]:
1507     raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
1508     return arg.trace(axis_offset)
1509 gross 429 elif isinstance(arg,float):
1510 gross 785 raise TypeError,"illegal argument type float."
1511 gross 429 elif isinstance(arg,int):
1512 gross 785 raise TypeError,"illegal argument type int."
1513 gross 429 else:
1514 gross 785 raise TypeError,"Unknown argument type."
1515 gross 429
1516 gross 492 def transpose(arg,axis_offset=None):
1517     """
1518 caltinay 3507 Returns the transpose of ``arg`` by swapping the first ``axis_offset`` and
1519     the last ``rank-axis_offset`` components.
1520 gross 492
1521 jfenwick 2625 :param arg: argument
1522 caltinay 3507 :type arg: `escript.Data`, `Symbol`, ``numpy.ndarray``, ``float``, ``int``
1523 jfenwick 2625 :param axis_offset: the first ``axis_offset`` components are swapped with the
1524     rest. ``axis_offset`` must be non-negative and less or
1525     equal to the rank of ``arg``. If ``axis_offset`` is not
1526     present ``int(r/2)`` where r is the rank of ``arg`` is
1527 caltinay 2169 used.
1528 jfenwick 2625 :type axis_offset: ``int``
1529     :return: transpose of ``arg``
1530 caltinay 3507 :rtype: `escript.Data`, `Symbol`, ``numpy.ndarray``, ``float``, ``int``
1531 jfenwick 2625 depending on the type of ``arg``
1532 gross 492 """
1533 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1534     if axis_offset==None: axis_offset=int(arg.ndim/2)
1535     return numpy.transpose(arg,axes=range(axis_offset,arg.ndim)+range(0,axis_offset))
1536 gross 492 elif isinstance(arg,escript.Data):
1537 gross 785 r=arg.getRank()
1538     if axis_offset==None: axis_offset=int(r/2)
1539     if axis_offset<0 or axis_offset>r:
1540     raise ValueError,"axis_offset must be between 0 and %s"%r
1541     return arg._transpose(axis_offset)
1542 gross 492 elif isinstance(arg,float):
1543 caltinay 2169 if not ( axis_offset==0 or axis_offset==None):
1544 gross 785 raise ValueError,"axis_offset must be 0 for float argument"
1545 gross 492 return arg
1546     elif isinstance(arg,int):
1547 caltinay 2169 if not ( axis_offset==0 or axis_offset==None):
1548 gross 785 raise ValueError,"axis_offset must be 0 for int argument"
1549 gross 492 return float(arg)
1550 caltinay 3507 elif isinstance(arg,Symbol):
1551     r=arg.getRank()
1552     if axis_offset==None: axis_offset=int(r/2)
1553     if axis_offset<0 or axis_offset>r:
1554     raise ValueError,"axis_offset must be between 0 and %s"%r
1555     return arg.transpose(axis_offset)
1556 gross 492 else:
1557 gross 785 raise TypeError,"Unknown argument type."
1558 gross 492
1559 gross 804 def swap_axes(arg,axis0=0,axis1=1):
1560     """
1561 jfenwick 2625 Returns the swap of ``arg`` by swapping the components ``axis0`` and ``axis1``.
1562 gross 804
1563 jfenwick 2625 :param arg: argument
1564 caltinay 3507 :type arg: `escript.Data`, `Symbol`, ``numpy.ndarray``
1565 jfenwick 2625 :param axis0: first axis. ``axis0`` must be non-negative and less than the
1566     rank of ``arg``.
1567     :type axis0: ``int``
1568     :param axis1: second axis. ``axis1`` must be non-negative and less than the
1569     rank of ``arg``.
1570     :type axis1: ``int``
1571     :return: ``arg`` with swapped components
1572 caltinay 3507 :rtype: `escript.Data`, `Symbol` or ``numpy.ndarray`` depending on the type
1573     of ``arg``
1574 gross 804 """
1575     if axis0 > axis1:
1576     axis0,axis1=axis1,axis0
1577 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1578     return numpy.swapaxes(arg,axis0,axis1)
1579 gross 804 elif isinstance(arg,escript.Data):
1580     return arg._swap_axes(axis0,axis1)
1581 caltinay 3507 elif isinstance(arg,Symbol):
1582     return arg.swap_axes(axis0,axis1)
1583 gross 804 elif isinstance(arg,float):
1584 caltinay 2169 raise TypeError,"float argument is not supported."
1585 gross 804 elif isinstance(arg,int):
1586 caltinay 2169 raise TypeError,"int argument is not supported."
1587 gross 804 else:
1588     raise TypeError,"Unknown argument type."
1589    
1590 gross 536 def symmetric(arg):
1591     """
1592 jfenwick 2625 Returns the symmetric part of the square matrix ``arg``. That is,
1593     *(arg+transpose(arg))/2*.
1594 gross 492
1595 jfenwick 2625 :param arg: input matrix. Must have rank 2 or 4 and be square.
1596 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1597 jfenwick 2625 :return: symmetric part of ``arg``
1598 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1599 gross 536 """
1600 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1601     if arg.ndim==2:
1602 gross 550 if not (arg.shape[0]==arg.shape[1]):
1603 gross 785 raise ValueError,"argument must be square."
1604 jfenwick 2455 elif arg.ndim==4:
1605 gross 550 if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
1606 gross 785 raise ValueError,"argument must be square."
1607 gross 550 else:
1608 gross 785 raise ValueError,"rank 2 or 4 is required."
1609 gross 550 return (arg+transpose(arg))/2
1610     elif isinstance(arg,escript.Data):
1611 gross 785 if arg.getRank()==2:
1612     if not (arg.getShape()[0]==arg.getShape()[1]):
1613     raise ValueError,"argument must be square."
1614     return arg._symmetric()
1615     elif arg.getRank()==4:
1616     if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
1617     raise ValueError,"argument must be square."
1618     return arg._symmetric()
1619     else:
1620     raise ValueError,"rank 2 or 4 is required."
1621 caltinay 3507 elif isinstance(arg, Symbol):
1622     if arg.getRank()==2:
1623     if arg.getShape()[0]!=arg.getShape()[1]:
1624     raise ValueError,"symmetric: argument must be square."
1625     elif arg.getRank()==4:
1626     if arg.getShape()[0]!=arg.getShape()[2] or arg.getShape()[1]!=arg.getShape()[3]:
1627     raise ValueError,"symmetric: argument must be square."
1628     else:
1629     raise ValueError,"symmetric: rank 2 or 4 is required."
1630     return (arg+transpose(arg))/2
1631 gross 550 elif isinstance(arg,float):
1632     return arg
1633     elif isinstance(arg,int):
1634     return float(arg)
1635     else:
1636     raise TypeError,"symmetric: Unknown argument type."
1637 gross 536
1638     def nonsymmetric(arg):
1639     """
1640 jfenwick 2625 Returns the non-symmetric part of the square matrix ``arg``. That is,
1641     *(arg-transpose(arg))/2*.
1642 gross 536
1643 jfenwick 2625 :param arg: input matrix. Must have rank 2 or 4 and be square.
1644 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1645 jfenwick 2625 :return: non-symmetric part of ``arg``
1646 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1647 gross 536 """
1648 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1649     if arg.ndim==2:
1650 gross 550 if not (arg.shape[0]==arg.shape[1]):
1651     raise ValueError,"nonsymmetric: argument must be square."
1652 jfenwick 2455 elif arg.ndim==4:
1653 gross 550 if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
1654     raise ValueError,"nonsymmetric: argument must be square."
1655     else:
1656     raise ValueError,"nonsymmetric: rank 2 or 4 is required."
1657     return (arg-transpose(arg))/2
1658     elif isinstance(arg,escript.Data):
1659 gross 785 if arg.getRank()==2:
1660     if not (arg.getShape()[0]==arg.getShape()[1]):
1661     raise ValueError,"argument must be square."
1662     return arg._nonsymmetric()
1663     elif arg.getRank()==4:
1664     if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
1665     raise ValueError,"argument must be square."
1666     return arg._nonsymmetric()
1667     else:
1668     raise ValueError,"rank 2 or 4 is required."
1669 caltinay 3507 elif isinstance(arg, Symbol):
1670     if arg.getRank()==2:
1671     if arg.getShape()[0]!=arg.getShape()[1]:
1672     raise ValueError,"nonsymmetric: argument must be square."
1673     elif arg.getRank()==4:
1674     if arg.getShape()[0]!=arg.getShape()[2] or arg.getShape()[1]!=arg.getShape()[3]:
1675     raise ValueError,"nonsymmetric: argument must be square."
1676     else:
1677     raise ValueError,"nonsymmetric: rank 2 or 4 is required."
1678     return (arg-transpose(arg))/2
1679 gross 550 elif isinstance(arg,float):
1680 caltinay 3507 return arg
1681 gross 550 elif isinstance(arg,int):
1682 caltinay 3507 return float(arg)
1683 gross 550 else:
1684 caltinay 3507 raise TypeError,"nonsymmetric: Unknown argument type."
1685 gross 536
1686 gross 433 def inverse(arg):
1687     """
1688 jfenwick 2625 Returns the inverse of the square matrix ``arg``.
1689 gross 433
1690 jfenwick 2625 :param arg: square matrix. Must have rank 2 and the first and second
1691 caltinay 2169 dimension must be equal.
1692 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1693 jfenwick 2625 :return: inverse of the argument. ``matrix_mult(inverse(arg),arg)`` will be
1694     almost equal to ``kronecker(arg.getShape()[0])``
1695 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1696 jfenwick 2625 :note: for `escript.Data` objects the dimension is restricted to 3.
1697 gross 433 """
1698 jfenwick 2455 import numpy.linalg
1699     if isinstance(arg,numpy.ndarray):
1700     return numpy.linalg.tensorinv(arg,ind=1)
1701 gross 433 elif isinstance(arg,escript.Data):
1702     return escript_inverse(arg)
1703     elif isinstance(arg,float):
1704     return 1./arg
1705     elif isinstance(arg,int):
1706     return 1./float(arg)
1707 caltinay 3507 elif isinstance(arg,Symbol):
1708 caltinay 3517 return arg.inverse()
1709 gross 433 else:
1710     raise TypeError,"inverse: Unknown argument type."
1711    
1712     def escript_inverse(arg): # this should be escript._inverse and use LAPACK
1713 caltinay 3493 "arg is a Data object!"
1714 jfenwick 2742 return arg._inverse()
1715 gross 433
1716 jfenwick 2742 #if not arg.getRank()==2:
1717     #raise ValueError,"escript_inverse: argument must have rank 2"
1718     #s=arg.getShape()
1719     #if not s[0] == s[1]:
1720     #raise ValueError,"escript_inverse: argument must be a square matrix."
1721     #out=escript.Data(0.,s,arg.getFunctionSpace())
1722     #if s[0]==1:
1723     #if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
1724     #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1725     #out[0,0]=1./arg[0,0]
1726     #elif s[0]==2:
1727     #A11=arg[0,0]
1728     #A12=arg[0,1]
1729     #A21=arg[1,0]
1730     #A22=arg[1,1]
1731     #D = A11*A22-A12*A21
1732     #if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
1733     #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1734     #D=1./D
1735     #out[0,0]= A22*D
1736     #out[1,0]=-A21*D
1737     #out[0,1]=-A12*D
1738     #out[1,1]= A11*D
1739     #elif s[0]==3:
1740     #A11=arg[0,0]
1741     #A21=arg[1,0]
1742     #A31=arg[2,0]
1743     #A12=arg[0,1]
1744     #A22=arg[1,1]
1745     #A32=arg[2,1]
1746     #A13=arg[0,2]
1747     #A23=arg[1,2]
1748     #A33=arg[2,2]
1749     #D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
1750     #if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
1751     #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1752     #D=1./D
1753     #out[0,0]=(A22*A33-A23*A32)*D
1754     #out[1,0]=(A31*A23-A21*A33)*D
1755     #out[2,0]=(A21*A32-A31*A22)*D
1756     #out[0,1]=(A13*A32-A12*A33)*D
1757     #out[1,1]=(A11*A33-A31*A13)*D
1758     #out[2,1]=(A12*A31-A11*A32)*D
1759     #out[0,2]=(A12*A23-A13*A22)*D
1760     #out[1,2]=(A13*A21-A11*A23)*D
1761     #out[2,2]=(A11*A22-A12*A21)*D
1762     #else:
1763     #raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
1764     #return out
1765    
1766 gross 528 def eigenvalues(arg):
1767     """
1768 jfenwick 2625 Returns the eigenvalues of the square matrix ``arg``.
1769 gross 528
1770 jfenwick 2625 :param arg: square matrix. Must have rank 2 and the first and second
1771 caltinay 2169 dimension must be equal. It must also be symmetric, ie.
1772 jfenwick 2625 ``transpose(arg)==arg`` (this is not checked).
1773 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1774 jfenwick 2625 :return: the eigenvalues in increasing order
1775 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1776     :note: for `escript.Data` and `Symbol` objects the dimension is
1777     restricted to 3.
1778 gross 528 """
1779 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1780     out=numpy.linalg.eigvals((arg+numpy.transpose(arg))/2.)
1781 gross 580 out.sort()
1782 gross 529 return out
1783 gross 530 elif isinstance(arg,escript.Data):
1784 gross 580 return arg._eigenvalues()
1785 gross 528 elif isinstance(arg,float):
1786 gross 529 return arg
1787 gross 528 elif isinstance(arg,int):
1788 gross 529 return float(arg)
1789 caltinay 3507 elif isinstance(arg,Symbol):
1790 caltinay 3496 return symfn.eigenvalues(arg)
1791 gross 528 else:
1792     raise TypeError,"eigenvalues: Unknown argument type."
1793 gross 587
1794     def eigenvalues_and_eigenvectors(arg):
1795     """
1796 jfenwick 2625 Returns the eigenvalues and eigenvectors of the square matrix ``arg``.
1797 gross 587
1798 jfenwick 2625 :param arg: square matrix. Must have rank 2 and the first and second
1799 caltinay 2169 dimension must be equal. It must also be symmetric, ie.
1800 jfenwick 2625 ``transpose(arg)==arg`` (this is not checked).
1801     :type arg: `escript.Data`
1802     :return: the eigenvalues and eigenvectors. The eigenvalues are ordered by
1803 caltinay 2169 increasing value. The eigenvectors are orthogonal and normalized.
1804     If V are the eigenvectors then V[:,i] is the eigenvector
1805     corresponding to the i-th eigenvalue.
1806 jfenwick 2625 :rtype: `tuple` of `escript.Data`
1807     :note: The dimension is restricted to 3.
1808 gross 587 """
1809 jfenwick 2455 if isinstance(arg,numpy.ndarray):
1810     raise TypeError,"eigenvalues_and_eigenvectors does not support numpy.ndarray arguments"
1811 gross 587 elif isinstance(arg,escript.Data):
1812     return arg._eigenvalues_and_eigenvectors()
1813     elif isinstance(arg,float):
1814 jfenwick 2455 return (numpy.array([[arg]],numpy.float_),numpy.ones((1,1),numpy.float_))
1815 gross 587 elif isinstance(arg,int):
1816 jfenwick 2455 return (numpy.array([[arg]],numpy.float_),numpy.ones((1,1),numpy.float_))
1817 caltinay 3507 elif isinstance(arg,Symbol):
1818 caltinay 3496 return symfn.eigenvalues_and_eigenvectors(arg)
1819 gross 587 else:
1820     raise TypeError,"eigenvalues: Unknown argument type."
1821 caltinay 2169
1822 caltinay 3609 def mult(arg0,arg1):
1823     """
1824     Product of ``arg0`` and ``arg1``.
1825    
1826     :param arg0: first term
1827     :type arg0: `Symbol`, ``float``, ``int``, `escript.Data` or
1828     ``numpy.ndarray``
1829     :param arg1: second term
1830     :type arg1: `Symbol`, ``float``, ``int``, `escript.Data` or
1831     ``numpy.ndarray``
1832     :return: the product of ``arg0`` and ``arg1``
1833     :rtype: `Symbol`, ``float``, ``int``, `escript.Data` or
1834     ``numpy.ndarray``
1835     :note: The shape of both arguments is matched according to the rules
1836     used in `matchShape`.
1837     """
1838     args=matchShape(arg0,arg1)
1839     if testForZero(args[0]) or testForZero(args[1]):
1840     return numpy.zeros(getShape(args[0]),numpy.float64)
1841     else:
1842     if isinstance(args[0],numpy.ndarray):
1843     return args[1]*args[0]
1844     else:
1845     return args[0]*args[1]
1846    
1847 gross 312 def maximum(*args):
1848 jgs 149 """
1849 jfenwick 2625 The maximum over arguments ``args``.
1850 caltinay 2169
1851 jfenwick 2625 :param args: arguments
1852 caltinay 3507 :type args: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1853     ``float``
1854 jfenwick 2625 :return: an object which in each entry gives the maximum of the
1855     corresponding values in ``args``
1856 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1857 jfenwick 2625 ``float`` depending on the input
1858 jgs 149 """
1859 caltinay 3507 if max([isinstance(v,Symbol) for v in args]):
1860 caltinay 3493 return symfn.maximum(*args)
1861 gross 312 out=None
1862     for a in args:
1863     if out==None:
1864 gross 2237 out=a*1.
1865 gross 312 else:
1866 gross 2237 if isinstance(out,escript.Data) and isinstance(a,escript.Data):
1867 caltinay 3493 if out.getRank()==0 and a.getRank()>0:
1868     # We need to consider the case where we have scalars and
1869     # higher ranked objects mixed. If the scalar was first it will
1870     # get picked as the initial out and we have a problem, so we
1871     # swap the objects
1872     res=a.copy() #Deep copy of a
1873     res.copyWithMask(out,wherePositive(out-a))
1874     out=res
1875     else:
1876     out.copyWithMask(a,wherePositive(a-out))
1877 gross 2237 else:
1878 caltinay 3609 if isinstance(a, numpy.ndarray):
1879 caltinay 3495 diff=-out+a
1880     else:
1881     diff=a-out
1882 caltinay 3609 temp=mult(whereNonPositive(diff),out)+mult(wherePositive(diff),a)
1883     if isinstance(out,numpy.ndarray) and isinstance(a,numpy.ndarray):
1884     # we need to convert the result to an array
1885     temp=numpy.array(temp)
1886     out=temp
1887 gross 312 return out
1888 caltinay 2169
1889 gross 396 def minimum(*args):
1890 jgs 149 """
1891 jfenwick 2625 The minimum over arguments ``args``.
1892 caltinay 2169
1893 jfenwick 2625 :param args: arguments
1894 caltinay 3507 :type args: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1895     ``float``
1896 jfenwick 2625 :return: an object which gives in each entry the minimum of the
1897     corresponding values in ``args``
1898 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1899     ``float`` depending on the input
1900 jgs 149 """
1901 caltinay 3507 if max([isinstance(v,Symbol) for v in args]):
1902 caltinay 3493 return symfn.minimum(*args)
1903 gross 312 out=None
1904     for a in args:
1905     if out==None:
1906 caltinay 3495 out=a*1.
1907 gross 312 else:
1908 gross 2237 if isinstance(out,escript.Data) and isinstance(a,escript.Data):
1909 caltinay 3493 if out.getRank()==0 and a.getRank()>0:
1910     # We need to consider the case where we have scalars and
1911     # higher ranked objects mixed. If the scalar was first it will
1912     # get picked as the initial out and we have a problem, so we
1913     # swap the objects
1914     res=a.copy() # Deep copy of a
1915     res.copyWithMask(out,whereNegative(out-a))
1916     out=res
1917     else:
1918     out.copyWithMask(a,whereNegative(a-out))
1919 gross 2237 else:
1920 caltinay 3609 if isinstance(a, numpy.ndarray):
1921 caltinay 3495 diff=-out+a
1922     else:
1923     diff=a-out
1924 caltinay 3609 #out=add(out,mult(whereNegative(diff),diff))
1925     temp=mult(whereNonNegative(diff),out)+mult(whereNegative(diff),a)
1926     if isinstance(out,numpy.ndarray) and isinstance(a,numpy.ndarray):
1927     # we need to convert the result to an array
1928     temp=numpy.array(temp)
1929     out=temp
1930 gross 312 return out
1931 gross 378
1932 gross 881 def clip(arg,minval=None,maxval=None):
1933 gross 378 """
1934 jfenwick 2625 Cuts the values of ``arg`` between ``minval`` and ``maxval``.
1935 caltinay 2169
1936 jfenwick 2625 :param arg: argument
1937 caltinay 3507 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1938     ``float``
1939 jfenwick 2625 :param minval: lower range. If None no lower range is applied
1940     :type minval: ``float`` or ``None``
1941     :param maxval: upper range. If None no upper range is applied
1942     :type maxval: ``float`` or ``None``
1943     :return: an object that contains all values from ``arg`` between ``minval``
1944     and ``maxval``
1945 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1946     ``float`` depending on the input
1947 jfenwick 2625 :raise ValueError: if ``minval>maxval``
1948 gross 378 """
1949 caltinay 3507 if isinstance(arg, Symbol):
1950     clip_item=lambda item: symfn.clip(item, minval, maxval)
1951     return arg.applyfunc(clip_item)
1952 gross 881 if not minval==None and not maxval==None:
1953     if minval>maxval:
1954 caltinay 2169 raise ValueError,"minval = %s must be less than maxval %s"%(minval,maxval)
1955 gross 881 if minval == None:
1956     tmp=arg
1957     else:
1958     tmp=maximum(minval,arg)
1959     if maxval == None:
1960     return tmp
1961     else:
1962     return minimum(tmp,maxval)
1963 gross 378
1964 gross 290 def inner(arg0,arg1):
1965     """
1966 caltinay 2169 Inner product of the two arguments. The inner product is defined as:
1967 jgs 123
1968 jfenwick 2625 C{out=Sigma_s arg0[s]*arg1[s]}
1969 gross 290
1970 jfenwick 2625 where s runs through ``arg0.Shape``.
1971 gross 290
1972 jfenwick 2625 ``arg0`` and ``arg1`` must have the same shape.
1973 caltinay 2169
1974 jfenwick 2625 :param arg0: first argument
1975 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
1976 jfenwick 2625 :param arg1: second argument
1977 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
1978 jfenwick 2625 :return: the inner product of ``arg0`` and ``arg1`` at each data point
1979 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``
1980 caltinay 2169 depending on the input
1981 jfenwick 2625 :raise ValueError: if the shapes of the arguments are not identical
1982 jgs 82 """
1983 gross 912 sh0=getShape(arg0)
1984     sh1=getShape(arg1)
1985 gross 290 if not sh0==sh1:
1986     raise ValueError,"inner: shape of arguments does not match"
1987 gross 429 return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
1988 jgs 82
1989 gross 785 def outer(arg0,arg1):
1990 jgs 82 """
1991 caltinay 2169 The outer product of the two arguments. The outer product is defined as:
1992 gross 785
1993 jfenwick 2625 ``out[t,s]=arg0[t]*arg1[s]``
1994 gross 785
1995 caltinay 2169 where
1996 jfenwick 2625 - s runs through ``arg0.Shape``
1997     - t runs through ``arg1.Shape``
1998 gross 785
1999 jfenwick 2625 :param arg0: first argument
2000 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2001 jfenwick 2625 :param arg1: second argument
2002 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2003 jfenwick 2625 :return: the outer product of ``arg0`` and ``arg1`` at each data point
2004 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2005 gross 785 """
2006     return generalTensorProduct(arg0,arg1,axis_offset=0)
2007    
2008 ksteube 795 def matrixmult(arg0,arg1):
2009 gross 785 """
2010 jfenwick 2625 See `matrix_mult`.
2011 gross 785 """
2012     return matrix_mult(arg0,arg1)
2013    
2014     def matrix_mult(arg0,arg1):
2015     """
2016 caltinay 2169 matrix-matrix or matrix-vector product of the two arguments.
2017 jgs 82
2018 jfenwick 2625 C{out[s0]=Sigma_{r0} arg0[s0,r0]*arg1[r0]}
2019 gross 290
2020 caltinay 2169 or
2021 gross 290
2022 jfenwick 2625 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[r0,s1]}
2023 caltinay 2169
2024 jfenwick 2625 The second dimension of ``arg0`` and the first dimension of ``arg1`` must
2025 caltinay 2169 match.
2026    
2027 jfenwick 2625 :param arg0: first argument of rank 2
2028 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2029 jfenwick 2625 :param arg1: second argument of at least rank 1
2030 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2031 jfenwick 2625 :return: the matrix-matrix or matrix-vector product of ``arg0`` and ``arg1``
2032 caltinay 2169 at each data point
2033 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2034 jfenwick 2625 :raise ValueError: if the shapes of the arguments are not appropriate
2035 jgs 147 """
2036 gross 912 sh0=getShape(arg0)
2037     sh1=getShape(arg1)
2038 gross 290 if not len(sh0)==2 :
2039     raise ValueError,"first argument must have rank 2"
2040     if not len(sh1)==2 and not len(sh1)==1:
2041     raise ValueError,"second argument must have rank 1 or 2"
2042 gross 429 return generalTensorProduct(arg0,arg1,axis_offset=1)
2043 jgs 147
2044 gross 785 def tensormult(arg0,arg1):
2045 jgs 147 """
2046 jfenwick 2625 See `tensor_mult`.
2047 jgs 149 """
2048 gross 785 return tensor_mult(arg0,arg1)
2049 jgs 149
2050 gross 785 def tensor_mult(arg0,arg1):
2051 jgs 149 """
2052 caltinay 2169 The tensor product of the two arguments.
2053 gross 290
2054 jfenwick 2625 For ``arg0`` of rank 2 this is
2055 gross 290
2056 jfenwick 2625 C{out[s0]=Sigma_{r0} arg0[s0,r0]*arg1[r0]}
2057 gross 290
2058 caltinay 2169 or
2059 jgs 149
2060 jfenwick 2625 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[r0,s1]}
2061 caltinay 2169
2062 jfenwick 2625 and for ``arg0`` of rank 4 this is
2063 caltinay 2169
2064 jfenwick 2625 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]}
2065 caltinay 2169
2066 gross 720 or
2067 gross 290
2068 jfenwick 2625 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]}
2069 gross 290
2070 caltinay 2169 or
2071 gross 290
2072 jfenwick 2625 C{out[s0,s1]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]}
2073 gross 290
2074 jfenwick 2625 In the first case the second dimension of ``arg0`` and the last dimension of
2075     ``arg1`` must match and in the second case the two last dimensions of ``arg0``
2076     must match the two first dimensions of ``arg1``.
2077 caltinay 2169
2078 jfenwick 2625 :param arg0: first argument of rank 2 or 4
2079 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2080 jfenwick 2625 :param arg1: second argument of shape greater than 1 or 2 depending on the
2081     rank of ``arg0``
2082 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2083 jfenwick 2625 :return: the tensor product of ``arg0`` and ``arg1`` at each data point
2084 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2085 gross 290 """
2086 gross 912 sh0=getShape(arg0)
2087     sh1=getShape(arg1)
2088 gross 290 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2089 gross 429 return generalTensorProduct(arg0,arg1,axis_offset=1)
2090 gross 290 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2091 gross 429 return generalTensorProduct(arg0,arg1,axis_offset=2)
2092 gross 290 else:
2093 caltinay 3493 raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
2094 gross 290
2095 gross 429 def generalTensorProduct(arg0,arg1,axis_offset=0):
2096 gross 290 """
2097 caltinay 2169 Generalized tensor product.
2098 gross 290
2099 jfenwick 2625 C{out[s,t]=Sigma_r arg0[s,r]*arg1[r,t]}
2100 gross 290
2101 gross 720 where
2102 jfenwick 2625 - s runs through ``arg0.Shape[:arg0.ndim-axis_offset]``
2103 caltinay 3609 - r runs through ``arg1.Shape[:axis_offset]``
2104 jfenwick 2625 - t runs through ``arg1.Shape[axis_offset:]``
2105 gross 290
2106 jfenwick 2625 :param arg0: first argument
2107 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2108 jfenwick 2625 :param arg1: second argument
2109 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2110 jfenwick 2625 :return: the general tensor product of ``arg0`` and ``arg1`` at each data
2111 caltinay 2169 point
2112 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2113 gross 290 """
2114     if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2115     arg0,arg1=matchType(arg0,arg1)
2116 caltinay 3507 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data,
2117     # or one is a Symbol and the other either of the allowed types
2118     if isinstance(arg0,Symbol):
2119     sh0=arg0.getShape()
2120     sh1=getShape(arg1)
2121     if not sh0[arg0.getRank()-axis_offset:]==sh1[:axis_offset]:
2122     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))
2123     if isinstance(arg1,float):
2124     return arg0*arg1
2125     elif isinstance(arg1,numpy.ndarray) or isinstance(arg1, Symbol):
2126 caltinay 3509 return arg0.tensorProduct(arg1, axis_offset)
2127 caltinay 3507 elif isinstance(arg1, escript.Data):
2128     raise TypeError("tensor product of Symbol and Data not supported yet")
2129     elif isinstance(arg0,numpy.ndarray):
2130 caltinay 3493 if not arg0.shape[arg0.ndim-axis_offset:]==arg1.shape[:axis_offset]:
2131     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)
2132     arg0_c=arg0.copy()
2133     arg1_c=arg1.copy()
2134     sh0,sh1=arg0.shape,arg1.shape
2135     d0,d1,d01=1,1,1
2136     for i in sh0[:arg0.ndim-axis_offset]: d0*=i
2137     for i in sh1[axis_offset:]: d1*=i
2138     for i in sh1[:axis_offset]: d01*=i
2139     arg0_c.resize((d0,d01))
2140     arg1_c.resize((d01,d1))
2141     out=numpy.zeros((d0,d1),numpy.float64)
2142     for i0 in range(d0):
2143     for i1 in range(d1):
2144     out[i0,i1]=numpy.sum(arg0_c[i0,:]*arg1_c[:,i1])
2145     out.resize(sh0[:arg0.ndim-axis_offset]+sh1[axis_offset:])
2146     return out
2147 gross 290 elif isinstance(arg0,escript.Data):
2148 caltinay 3509 if isinstance(arg1, Symbol):
2149     raise TypeError("tensor product of Data and Symbol not supported yet")
2150 caltinay 3493 return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2151 caltinay 3507 raise TypeError("generalTensorProduct: Unsupported argument type")
2152 caltinay 2169
2153 ksteube 813 def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
2154 caltinay 2169 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2155 ksteube 813 return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
2156 ksteube 795
2157 gross 785 def transposed_matrix_mult(arg0,arg1):
2158     """
2159 caltinay 2169 transposed(matrix)-matrix or transposed(matrix)-vector product of the two
2160     arguments.
2161 gross 785
2162 jfenwick 2625 C{out[s0]=Sigma_{r0} arg0[r0,s0]*arg1[r0]}
2163 gross 785
2164 caltinay 2169 or
2165 gross 785
2166 jfenwick 2625 C{out[s0,s1]=Sigma_{r0} arg0[r0,s0]*arg1[r0,s1]}
2167 gross 785
2168 jfenwick 2625 The function call ``transposed_matrix_mult(arg0,arg1)`` is equivalent to
2169     ``matrix_mult(transpose(arg0),arg1)``.
2170 caltinay 2169
2171 jfenwick 2625 The first dimension of ``arg0`` and ``arg1`` must match.
2172 caltinay 2169
2173 jfenwick 2625 :param arg0: first argument of rank 2
2174 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2175 jfenwick 2625 :param arg1: second argument of at least rank 1
2176 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2177 jfenwick 2625 :return: the product of the transpose of ``arg0`` and ``arg1`` at each data
2178 caltinay 2169 point
2179 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2180 jfenwick 2625 :raise ValueError: if the shapes of the arguments are not appropriate
2181 gross 785 """
2182 gross 912 sh0=getShape(arg0)
2183     sh1=getShape(arg1)
2184 gross 785 if not len(sh0)==2 :
2185     raise ValueError,"first argument must have rank 2"
2186     if not len(sh1)==2 and not len(sh1)==1:
2187     raise ValueError,"second argument must have rank 1 or 2"
2188     return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
2189    
2190     def transposed_tensor_mult(arg0,arg1):
2191     """
2192 caltinay 2169 The tensor product of the transpose of the first and the second argument.
2193 gross 785
2194 jfenwick 2625 For ``arg0`` of rank 2 this is
2195 gross 785
2196 jfenwick 2625 C{out[s0]=Sigma_{r0} arg0[r0,s0]*arg1[r0]}
2197 gross 785
2198 caltinay 2169 or
2199 gross 785
2200 jfenwick 2625 C{out[s0,s1]=Sigma_{r0} arg0[r0,s0]*arg1[r0,s1]}
2201 caltinay 2169
2202 jfenwick 2625 and for ``arg0`` of rank 4 this is
2203 caltinay 2169
2204 jfenwick 2625 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]}
2205 caltinay 2169
2206 gross 785 or
2207    
2208 jfenwick 2625 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]}
2209 gross 785
2210 caltinay 2169 or
2211 gross 785
2212 jfenwick 2625 C{out[s0,s1]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]}
2213 gross 785
2214 jfenwick 2625 In the first case the first dimension of ``arg0`` and the first dimension of
2215     ``arg1`` must match and in the second case the two first dimensions of
2216     ``arg0`` must match the two first dimensions of ``arg1``.
2217 gross 785
2218 jfenwick 2625 The function call ``transposed_tensor_mult(arg0,arg1)`` is equivalent to
2219     ``tensor_mult(transpose(arg0),arg1)``.
2220 caltinay 2169
2221 jfenwick 2625 :param arg0: first argument of rank 2 or 4
2222 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2223 jfenwick 2625 :param arg1: second argument of shape greater of 1 or 2 depending on the
2224     rank of ``arg0``
2225 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2226 jfenwick 2625 :return: the tensor product of transpose of arg0 and arg1 at each data point
2227 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2228 gross 785 """
2229 gross 912 sh0=getShape(arg0)
2230     sh1=getShape(arg1)
2231 gross 785 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2232     return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
2233     elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2234     return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
2235     else:
2236     raise ValueError,"first argument must have rank 2 or 4"
2237    
2238     def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
2239     """
2240 jfenwick 2625 Generalized tensor product of transposed of ``arg0`` and ``arg1``.
2241 gross 785
2242 jfenwick 2625 C{out[s,t]=Sigma_r arg0[r,s]*arg1[r,t]}
2243 gross 785
2244     where
2245 jfenwick 2625 - s runs through ``arg0.Shape[axis_offset:]``
2246     - r runs through ``arg0.Shape[:axis_offset]``
2247     - t runs through ``arg1.Shape[axis_offset:]``
2248 gross 785
2249 jfenwick 2625 The function call ``generalTransposedTensorProduct(arg0,arg1,axis_offset)``
2250 caltinay 2169 is equivalent to
2251 jfenwick 2625 ``generalTensorProduct(transpose(arg0,arg0.ndim-axis_offset),arg1,axis_offset)``.
2252 gross 785
2253 jfenwick 2625 :param arg0: first argument
2254 caltinay 3507 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2255 jfenwick 2625 :param arg1: second argument
2256 caltinay 3507 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2257 jfenwick 2625 :return: the general tensor product of ``transpose(arg0)`` and ``arg1`` at
2258 caltinay 2169 each data point
2259 caltinay 3507 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2260 gross 785 """
2261     if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2262     arg0,arg1=matchType(arg0,arg1)
2263 caltinay 3509 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data,
2264     # or one is a Symbol and the other either of the allowed types
2265     if isinstance(arg0,Symbol):
2266     sh0=arg0.getShape()
2267     sh1=getShape(arg1)
2268     if not sh0[:axis_offset]==sh1[:axis_offset]:
2269     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))
2270     if isinstance(arg1,float):
2271     return arg0*arg1
2272     elif isinstance(arg1,numpy.ndarray) or isinstance(arg1, Symbol):
2273     return arg0.transposedTensorProduct(arg1, axis_offset)
2274     elif isinstance(arg1, escript.Data):
2275     raise TypeError("tensor product of Symbol and Data not supported yet")
2276     elif isinstance(arg0,numpy.ndarray):
2277 caltinay 3493 if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
2278     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)
2279     arg0_c=arg0.copy()
2280     arg1_c=arg1.copy()
2281     sh0,sh1=arg0.shape,arg1.shape
2282     d0,d1,d01=1,1,1
2283     for i in sh0[axis_offset:]: d0*=i
2284     for i in sh1[axis_offset:]: d1*=i
2285     for i in sh0[:axis_offset]: d01*=i
2286     arg0_c.resize((d01,d0))
2287     arg1_c.resize((d01,d1))
2288     out=numpy.zeros((d0,d1),numpy.float64)
2289     for i0 in range(d0):
2290     for i1 in range