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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 105323 byte(s)
Fast forward to latest trunk revision 3788.

1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # Primary Business: Queensland, Australia
10 # Licensed under the Open Software License version 3.0
11 # http://www.opensource.org/licenses/osl-3.0.php
12 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 Utility functions for escript
25
26 :var __author__: name of author
27 :var __copyright__: copyrights
28 :var __license__: licence agreement
29 :var __url__: url entry point on documentation
30 :var __version__: version
31 :var __date__: date of the version
32 :var EPSILON: smallest positive value with 1.<1.+EPSILON
33 :var DBLE_MAX: largest positive float
34 """
35
36 __author__="Lutz Gross, l.gross@uq.edu.au"
37
38
39 import math
40 import numpy
41 from . import escript
42 import os
43 from esys.escript import C_GeneralTensorProduct
44 from esys.escript import getVersion, getMPIRankWorld, getMPIWorldMax
45 from esys.escript import printParallelThreadCounts
46 from esys.escript import listEscriptParams
47 from esys.escript.escriptcpp import Data, _saveDataCSV, _condEval
48 from esys.escript.symbolic import *
49
50 #=========================================================
51 # some helpers:
52 #=========================================================
53 def getEpsilon():
54 return escript.getMachinePrecision()
55 EPSILON=getEpsilon()
56
57 def getMaxFloat():
58 return escript.getMaxFloat()
59 DBLE_MAX=getMaxFloat()
60
61 def getTagNames(domain):
62 """
63 Returns a list of tag names used by the domain.
64
65 :param domain: a domain object
66 :type domain: `escript.Domain`
67 :return: a list of tag names used by the domain
68 :rtype: ``list`` of ``str``
69 """
70 return [n.strip() for n in domain.showTagNames().split(",") ]
71
72 def insertTagNames(domain,**kwargs):
73 """
74 Inserts tag names into the domain.
75
76 :param domain: a domain object
77 :type domain: ``escript.Domain``
78 :keyword <tag_name>: tag key assigned to <tag_name>
79 :type <tag_name>: ``int``
80 """
81 for k in kwargs:
82 domain.setTagMap(k,kwargs[k])
83
84 def insertTaggedValues(target,**kwargs):
85 """
86 Inserts tagged values into the target using tag names.
87
88 :param target: data to be filled by tagged values
89 :type target: `escript.Data`
90 :keyword <tag_name>: value to be used for <tag_name>
91 :type <tag_name>: ``float`` or ``numpy.ndarray``
92 :return: ``target``
93 :rtype: `escript.Data`
94 """
95 for k in kwargs:
96 target.setTaggedValue(k,kwargs[k])
97 return target
98
99
100 def interpolateTable(tab, dat, start, step, undef=1.e50, check_boundaries=False):
101 try:
102 dim=len(start)
103 except TypeError:
104 start=(start,)
105 dim=1
106 try:
107 slen=len(step)
108 except TypeError:
109 step=(step,)
110 slen=1
111 if dim<1 or dim>3:
112 raise ValueError("Length of start list must be between 1 and 3.")
113 if dim!=slen:
114 raise ValueError("Length of start and step must be the same.")
115 dshape=dat.getShape()
116 if len(dshape)==0:
117 datdim=0
118 firstdim=dat
119 else:
120 datdim=dshape[0]
121 firstdim=dat[0]
122 #So now we know firstdim is a scalar
123 if (dim==1 and datdim>1) or (dim>1 and datdim!=dim):
124 print((dim, datdim))
125 raise ValueError("The dimension of dat must be equal to the length of start.")
126 if dim==3:
127 d1=dat[1]
128 d2=dat[2]
129 return firstdim._interpolateTable3d(tab, start[0], step[0], d1, start[1], step[1], d2, start[2], step[2], undef, check_boundaries)
130 if dim==2:
131 d1=dat[1]
132 return firstdim.interpolateTable(tab, start[0], step[0], d1, start[1], step[1], undef, check_boundaries)
133 # return d1.interpolateTable(tab, start[1], step[1], firstdim, start[0], step[0], undef, check_boundaries)
134 else:
135 return firstdim.interpolateTable(tab, start[0], step[0], undef, check_boundaries)
136
137
138 def saveDataCSV(filename, append=False, sep=", ", csep="_", **data):
139 """
140 Writes `Data` objects to a csv file.
141 These objects must have compatible FunctionSpaces. ie it must be possible to
142 interpolate all data to one `FunctionSpace`.
143
144 :param filename: file to save data to.
145 :type filename: ``string``
146 :param append: If ``True``, then open file at end rather than beginning
147 :type append: ``bool``
148 :param sep: separator between fields
149 :type sep: ``string``
150 :param csep: separator for components of rank2 and above eg ('_' -> c0_1)
151
152 The keyword args are Data objects to save.
153 If a scalar `Data` object is passed with the name ``mask``, then only
154 samples which correspond to positive values in ``mask`` will be output.
155 Example::
156
157 s=Scalar(..)
158 v=Vector(..)
159 t=Tensor(..)
160 f=float()
161 saveDataCSV("f.csv",a=s, b=v, c=t, d=f)
162
163 Will result in a file
164
165 a, b0, b1, c0_0, c0_1, .., c1_0 d
166 1.0, 1.5, 2.7, 3.1, 3.4, .., 0.89 0.0
167 0.9, 8.7, 1.9, 3.4, 7.8, .., 1.21 0.0
168
169 The first line is a header, the remaining lines give the values.
170
171 """
172 # find a function space:
173 fs = None
174 for n,d in list(data.items()):
175 if isinstance(d, Data): fs=d.getFunctionSpace()
176 if fs == None:
177 raise ValueError("saveDataCSV: there must be at least one Data object in the argument list.")
178
179 new_data={}
180 for n,d in list(data.items()):
181 if isinstance(d, Data):
182 new_data[n]=d
183 else:
184 try:
185 new_data[n]=Data(d,fs)
186 except:
187 raise ValueError("saveDataCSV: unknown non-data argument type for %s"%(str(n)))
188 _saveDataCSV(filename, new_data,sep, csep, append)
189
190 def saveVTK(filename,domain=None, metadata=None, metadata_schema=None, **data):
191 """
192 Deprecated. See esys.weipa.saveVTK().
193 """
194
195 import warnings
196 msg = "esys.escript.util.saveVTK is deprecated. Use the esys.weipa module."
197 warnings.warn(msg, DeprecationWarning, stacklevel=2)
198
199 from esys import weipa
200 weipa.saveVTK(filename, domain, metadata, metadata_schema, **data)
201
202 def saveDX(filename,domain=None,**data):
203 """
204 Writes `Data` objects into a file using the OpenDX file format.
205
206 Example::
207
208 tmp=Scalar(..)
209 v=Vector(..)
210 saveDX("solution.dx", temperature=tmp, velocity=v)
211
212 ``tmp`` and ``v`` are written into "solution.dx" where ``tmp`` is named
213 "temperature" and ``v`` is named "velocity".
214
215 :param filename: file name of the output file
216 :type filename: ``str``
217 :param domain: domain of the `Data` objects. If not specified, the domain
218 of the given `Data` objects is used.
219 :type domain: `escript.Domain`
220 :keyword <name>: writes the assigned value to the DX file using <name> as
221 identifier. The identifier can be used to select the data
222 set when data are imported into DX.
223 :type <name>: `Data` object
224 :note: The data objects have to be defined on the same domain. They may not
225 be in the same `FunctionSpace` but one cannot expect that all
226 `FunctionSpace` s can be mixed. Typically, data on the boundary and
227 data on the interior cannot be mixed.
228 """
229 new_data={}
230 for n,d in list(data.items()):
231 if not d.isEmpty():
232 fs=d.getFunctionSpace()
233 domain2=fs.getDomain()
234 if fs == escript.Solution(domain2):
235 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
236 elif fs == escript.ReducedSolution(domain2):
237 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
238 elif fs == escript.ContinuousFunction(domain2):
239 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
240 else:
241 new_data[n]=d
242 if domain==None: domain=domain2
243 if domain==None:
244 raise ValueError("saveDX: no domain detected.")
245 domain.saveDX(filename,new_data)
246
247 def saveESD(datasetName, dataDir=".", domain=None, timeStep=0, deltaT=1, dynamicMesh=0, **data):
248 """
249 Saves `Data` objects to files and creates an I{escript dataset} (ESD) file
250 for convenient processing/visualisation.
251
252 Single timestep example::
253
254 tmp = Scalar(..)
255 v = Vector(..)
256 saveESD("solution", "data", temperature=tmp, velocity=v)
257
258 Time series example::
259
260 while t < t_end:
261 tmp = Scalar(..)
262 v = Vector(..)
263 # save every 10 timesteps
264 if t % 10 == 0:
265 saveESD("solution", "data", timeStep=t, deltaT=10, temperature=tmp, velocity=v)
266 t = t + 1
267
268 tmp, v and the domain are saved in native format in the "data"
269 directory and the file "solution.esd" is created that refers to tmp by
270 the name "temperature" and to v by the name "velocity".
271
272 :param datasetName: name of the dataset, used to name the ESD file
273 :type datasetName: ``str``
274 :param dataDir: optional directory where the data files should be saved
275 :type dataDir: ``str``
276 :param domain: domain of the `Data` object(s). If not specified, the
277 domain of the given `Data` objects is used.
278 :type domain: `escript.Domain`
279 :param timeStep: current timestep or sequence number - first one must be 0
280 :type timeStep: `int`
281 :param deltaT: timestep or sequence increment, see example above
282 :type deltaT: `int`
283 :param dynamicMesh: by default the mesh is assumed to be static and thus
284 only saved once at timestep 0 to save disk space.
285 Setting this to 1 changes the behaviour and the mesh
286 is saved at each timestep.
287 :type dynamicMesh: `int`
288 :keyword <name>: writes the assigned value to the file using <name> as
289 identifier
290 :type <name>: `Data` object
291 :note: The ESD concept is experimental and the file format likely to
292 change so use this function with caution.
293 :note: The data objects have to be defined on the same domain (but not
294 necessarily on the same `FunctionSpace`).
295 :note: When saving a time series the first timestep must be 0 and it is
296 assumed that data from all timesteps share the domain. The dataset
297 file is updated in each iteration.
298 """
299 new_data = {}
300 for n,d in list(data.items()):
301 if not d.isEmpty():
302 fs = d.getFunctionSpace()
303 domain2 = fs.getDomain()
304 if fs == escript.Solution(domain2):
305 new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
306 elif fs == escript.ReducedSolution(domain2):
307 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
308 else:
309 new_data[n]=d
310 if domain==None: domain=domain2
311 if domain==None:
312 raise ValueError("saveESD: no domain detected.")
313
314 if domain.onMasterProcessor() and not os.path.isdir(dataDir):
315 os.mkdir(dataDir)
316
317 meshFile = datasetName+"_mesh"
318 fileNumber = timeStep / deltaT
319
320 if dynamicMesh == 0:
321 # later timesteps reuse mesh from t=0
322 if timeStep == 0:
323 domain.dump(os.path.join(dataDir, meshFile + ".nc"))
324 else:
325 meshFile += ".%04d"
326 domain.dump(os.path.join(dataDir, (meshFile + ".nc") % fileNumber))
327
328 outputString = ""
329
330 if domain.onMasterProcessor():
331 outputString += "#escript datafile V1.0\n"
332 # number of timesteps (currently only 1 is supported)
333 outputString += "T=%d\n" % (fileNumber+1)
334 # timestep increment
335 outputString += "DT=1\n"
336 # name of the mesh file
337 outputString += "M=%s\n" % meshFile
338 # number of blocks (MPI size)
339 outputString += "N=%d\n" % domain.getMPISize()
340
341 # now add the variables
342 for varName, d in list(new_data.items()):
343 varFile = datasetName+"_"+varName+".%04d"
344 d.dump(os.path.join(dataDir, (varFile + ".nc") % fileNumber))
345 if domain.onMasterProcessor():
346 outputString += "V=%s:%s\n" % (varFile, varName)
347
348 if domain.onMasterProcessor():
349 esdfile = open(os.path.join(dataDir, datasetName+".esd"), "w")
350 esdfile.write(outputString)
351 esdfile.close()
352
353 def kronecker(d=3):
354 """
355 Returns the kronecker delta-symbol.
356
357 :param d: dimension or an object that has the ``getDim`` method defining the
358 dimension
359 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
360 :return: the object u of rank 2 with *u[i,j]=1* for *i=j* and *u[i,j]=0*
361 otherwise
362 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 2
363 """
364 return identityTensor(d)
365
366 def identity(shape=()):
367 """
368 Returns the ``shape`` x ``shape`` identity tensor.
369
370 :param shape: input shape for the identity tensor
371 :type shape: ``tuple`` of ``int``
372 :return: array whose shape is shape x shape where *u[i,k]=1* for *i=k* and
373 *u[i,k]=0* otherwise for len(shape)=1. If len(shape)=2:
374 *u[i,j,k,l]=1* for *i=k and j=l* and *u[i,j,k,l]=0* otherwise.
375 :rtype: ``numpy.ndarray`` of rank 1, rank 2 or rank 4
376 :raise ValueError: if len(shape)>2
377 """
378 if len(shape)>0:
379 out=numpy.zeros(shape+shape,numpy.float64)
380 if len(shape)==1:
381 for i0 in range(shape[0]):
382 out[i0,i0]=1.
383 elif len(shape)==2:
384 for i0 in range(shape[0]):
385 for i1 in range(shape[1]):
386 out[i0,i1,i0,i1]=1.
387 else:
388 raise ValueError("identity: length of shape is restricted to 2.")
389 else:
390 out=1.
391 return out
392
393 def zeros(shape=()):
394 """
395 Returns the ``shape`` zero tensor.
396
397 :param shape: input shape for the identity tensor
398 :type shape: ``tuple`` of ``int``
399 :return: array of shape filled with zeros
400 :rtype: ``numpy.ndarray``
401 """
402 if len(shape)>0:
403 out=numpy.zeros(shape,numpy.float64)
404 else:
405 out=0.
406 return out
407
408 def identityTensor(d=3):
409 """
410 Returns the ``d`` x ``d`` identity matrix.
411
412 :param d: dimension or an object that has the ``getDim`` method defining the
413 dimension
414 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
415 :return: the object u of rank 2 with *u[i,j]=1* for *i=j* and *u[i,j]=0*
416 otherwise
417 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 2
418 """
419 if isinstance(d,escript.FunctionSpace):
420 return escript.Data(identity((d.getDim(),)),d)
421 elif isinstance(d,escript.Domain):
422 return identity((d.getDim(),))
423 else:
424 return identity((d,))
425
426 def identityTensor4(d=3):
427 """
428 Returns the ``d`` x ``d`` x ``d`` x ``d`` identity tensor.
429
430 :param d: dimension or an object that has the ``getDim`` method defining the
431 dimension
432 :type d: ``int`` or any object with a ``getDim`` method
433 :return: the object u of rank 4 with *u[i,j,k,l]=1* for *i=k and j=l* and
434 *u[i,j,k,l]=0* otherwise
435 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 4
436 """
437 if isinstance(d,escript.FunctionSpace):
438 return escript.Data(identity((d.getDim(),d.getDim())),d)
439 elif isinstance(d,escript.Domain):
440 return identity((d.getDim(),d.getDim()))
441 else:
442 return identity((d,d))
443
444 def unitVector(i=0,d=3):
445 """
446 Returns a unit vector u of dimension d whose non-zero element is at index i.
447
448 :param i: index for non-zero element
449 :type i: ``int``
450 :param d: dimension or an object that has the ``getDim`` method defining the
451 dimension
452 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
453 :return: the object u of rank 1 with *u[j]=1* for *j=index* and *u[j]=0*
454 otherwise
455 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 1
456 """
457 return kronecker(d)[i]
458
459 #=========================================================================
460 # global reduction operations
461 #=========================================================================
462 def Lsup(arg):
463 """
464 Returns the Lsup-norm of argument ``arg``. This is the maximum absolute value
465 over all data points. This function is equivalent to ``sup(abs(arg))``.
466
467 :param arg: argument
468 :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
469 :return: maximum value of the absolute value of ``arg`` over all components
470 and all data points
471 :rtype: ``float``
472 :raise TypeError: if type of ``arg`` cannot be processed
473 """
474 if isinstance(arg,numpy.ndarray):
475 return sup(abs(arg))
476 elif isinstance(arg,escript.Data):
477 return arg._Lsup()
478 elif isinstance(arg,float):
479 return abs(arg)
480 elif isinstance(arg,int):
481 return abs(float(arg))
482 else:
483 raise TypeError("Lsup: Unknown argument type.")
484
485 def sup(arg):
486 """
487 Returns the maximum value over all data points.
488
489 :param arg: argument
490 :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
491 :return: maximum value of ``arg`` over all components and all data points
492 :rtype: ``float``
493 :raise TypeError: if type of ``arg`` cannot be processed
494 """
495 if isinstance(arg,numpy.ndarray):
496 return arg.max()
497 elif isinstance(arg,escript.Data):
498 return arg._sup()
499 elif isinstance(arg,float):
500 return arg
501 elif isinstance(arg,int):
502 return float(arg)
503 else:
504 raise TypeError("sup: Unknown argument type.")
505
506 def inf(arg):
507 """
508 Returns the minimum value over all data points.
509
510 :param arg: argument
511 :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
512 :return: minimum value of ``arg`` over all components and all data points
513 :rtype: ``float``
514 :raise TypeError: if type of ``arg`` cannot be processed
515 """
516 if isinstance(arg,numpy.ndarray):
517 return arg.min()
518 elif isinstance(arg,escript.Data):
519 return arg._inf()
520 elif isinstance(arg,float):
521 return arg
522 elif isinstance(arg,int):
523 return float(arg)
524 else:
525 raise TypeError("inf: Unknown argument type.")
526
527
528 #=========================================================================
529 # some little helpers
530 #=========================================================================
531 def getRank(arg):
532 """
533 Identifies the rank of the argument.
534
535 :param arg: an object whose rank is to be returned
536 :type arg: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``,
537 ``Symbol``
538 :return: the rank of the argument
539 :rtype: ``int``
540 :raise TypeError: if type of ``arg`` cannot be processed
541 """
542
543 if isinstance(arg,numpy.ndarray):
544 return arg.ndim
545 elif isinstance(arg,escript.Data):
546 return arg.getRank()
547 elif isinstance(arg,float):
548 return 0
549 elif isinstance(arg,int):
550 return 0
551 elif isinstance(arg,Symbol):
552 return arg.getRank()
553 else:
554 raise TypeError("getRank: Unknown argument type.")
555
556 def getShape(arg):
557 """
558 Identifies the shape of the argument.
559
560 :param arg: an object whose shape is to be returned
561 :type arg: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``,
562 ``Symbol``
563 :return: the shape of the argument
564 :rtype: ``tuple`` of ``int``
565 :raise TypeError: if type of ``arg`` cannot be processed
566 """
567
568 if isinstance(arg,numpy.ndarray):
569 return arg.shape
570 elif isinstance(arg,list):
571 return numpy.array(arg).shape
572 elif isinstance(arg,escript.Data):
573 return arg.getShape()
574 elif isinstance(arg,float):
575 return ()
576 elif isinstance(arg,int):
577 return ()
578 elif isinstance(arg,Symbol):
579 return arg.getShape()
580 else:
581 raise TypeError("getShape: Cannot identify shape")
582
583 def pokeDim(arg):
584 """
585 Identifies the spatial dimension of the argument.
586
587 :param arg: an object whose spatial dimension is to be returned
588 :type arg: any
589 :return: the spatial dimension of the argument, if available, or ``None``
590 :rtype: ``int`` or ``None``
591 """
592
593 if isinstance(arg,escript.Data):
594 return arg.getFunctionSpace().getDim()
595 else:
596 return None
597
598 def commonShape(arg0, arg1):
599 """
600 Returns a shape to which ``arg0`` can be extended from the right and ``arg1``
601 can be extended from the left.
602
603 :param arg0: an object with a shape (see `getShape`)
604 :param arg1: an object with a shape (see `getShape`)
605 :return: the shape of ``arg0`` or ``arg1`` such that the left part equals the
606 shape of ``arg0`` and the right end equals the shape of ``arg1``
607 :rtype: ``tuple`` of ``int``
608 :raise ValueError: if no shape can be found
609 """
610 sh0=getShape(arg0)
611 sh1=getShape(arg1)
612 if len(sh0)<len(sh1):
613 if not sh0==sh1[:len(sh0)]:
614 raise ValueError("argument 0 cannot be extended to the shape of argument 1")
615 return sh1
616 elif len(sh0)>len(sh1):
617 if not sh1==sh0[:len(sh1)]:
618 raise ValueError("argument 1 cannot be extended to the shape of argument 0")
619 return sh0
620 else:
621 if not sh0==sh1:
622 raise ValueError("argument 1 and argument 0 have not the same shape.")
623 return sh0
624
625 def commonDim(*args):
626 """
627 Identifies, if possible, the spatial dimension across a set of objects
628 which may or may not have a spatial dimension.
629
630 :param args: given objects
631 :return: the spatial dimension of the objects with identifiable dimension
632 (see `pokeDim`). If none of the objects has a spatial dimension
633 ``None`` is returned.
634 :rtype: ``int`` or ``None``
635 :raise ValueError: if the objects with identifiable dimension don't have
636 the same spatial dimension.
637 """
638 out=None
639 for a in args:
640 d=pokeDim(a)
641 if not out==None:
642 if not (d==None or out==d):
643 raise ValueError("dimension of arguments don't match")
644 else:
645 out=d
646 return out
647
648 def testForZero(arg):
649 """
650 Tests if the argument is identical to zero.
651
652 :param arg: the object to test for zero
653 :type arg: typically ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
654 :return: True if the argument is identical to zero, False otherwise
655 :rtype: ``bool``
656 """
657 if isinstance(arg,numpy.ndarray):
658 return not Lsup(arg)>0.
659 elif isinstance(arg,escript.Data):
660 return False
661 elif isinstance(arg,float):
662 return not Lsup(arg)>0.
663 elif isinstance(arg,int):
664 return not Lsup(arg)>0.
665 else:
666 return False
667
668 def matchType(arg0=0.,arg1=0.):
669 """
670 Converts ``arg0`` and ``arg1`` both to the same type ``numpy.ndarray`` or
671 `escript.Data`
672
673 :param arg0: first argument
674 :type arg0: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, ``Symbol``
675 :param arg1: second argument
676 :type arg1: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, ``Symbol``
677 :return: a tuple representing ``arg0`` and ``arg1`` with the same type or
678 with at least one of them being a `Symbol`
679 :rtype: ``tuple`` of two ``numpy.ndarray`` or two `escript.Data`
680 :raise TypeError: if type of ``arg0`` or ``arg1`` cannot be processed
681 """
682 if isinstance(arg0,numpy.ndarray):
683 if isinstance(arg1,numpy.ndarray):
684 pass
685 elif isinstance(arg1,escript.Data):
686 arg0=escript.Data(arg0,arg1.getFunctionSpace())
687 elif isinstance(arg1,float):
688 arg1=numpy.array(arg1,dtype=numpy.float64)
689 elif isinstance(arg1,int):
690 arg1=numpy.array(float(arg1),dtype=numpy.float64)
691 elif isinstance(arg1,Symbol):
692 pass
693 else:
694 raise TypeError("function: Unknown type of second argument.")
695 elif isinstance(arg0,escript.Data):
696 if isinstance(arg1,numpy.ndarray):
697 arg1=escript.Data(arg1,arg0.getFunctionSpace())
698 elif isinstance(arg1,escript.Data):
699 pass
700 elif isinstance(arg1,float):
701 arg1=escript.Data(arg1,(),arg0.getFunctionSpace())
702 elif isinstance(arg1,int):
703 arg1=escript.Data(float(arg1),(),arg0.getFunctionSpace())
704 elif isinstance(arg1,Symbol):
705 pass
706 else:
707 raise TypeError("function: Unknown type of second argument.")
708 elif isinstance(arg0,Symbol):
709 if isinstance(arg1,numpy.ndarray):
710 pass
711 elif isinstance(arg1,escript.Data):
712 pass
713 elif isinstance(arg1,float):
714 pass
715 elif isinstance(arg1,int):
716 pass
717 elif isinstance(arg1,Symbol):
718 pass
719 else:
720 raise TypeError("function: Unknown type of second argument.")
721 elif isinstance(arg0,float):
722 if isinstance(arg1,numpy.ndarray):
723 arg0=numpy.array(arg0,dtype=numpy.float64)
724 elif isinstance(arg1,escript.Data):
725 arg0=escript.Data(arg0,arg1.getFunctionSpace())
726 elif isinstance(arg1,float):
727 arg0=numpy.array(arg0,dtype=numpy.float64)
728 arg1=numpy.array(arg1,dtype=numpy.float64)
729 elif isinstance(arg1,int):
730 arg0=numpy.array(arg0,dtype=numpy.float64)
731 arg1=numpy.array(float(arg1),dtype=numpy.float64)
732 elif isinstance(arg1,Symbol):
733 pass
734 else:
735 raise TypeError("function: Unknown type of second argument.")
736 elif isinstance(arg0,int):
737 if isinstance(arg1,numpy.ndarray):
738 arg0=numpy.array(float(arg0),dtype=numpy.float64)
739 elif isinstance(arg1,escript.Data):
740 arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
741 elif isinstance(arg1,float):
742 arg0=numpy.array(float(arg0),dtype=numpy.float64)
743 arg1=numpy.array(arg1,dtype=numpy.float64)
744 elif isinstance(arg1,int):
745 arg0=numpy.array(float(arg0),dtype=numpy.float64)
746 arg1=numpy.array(float(arg1),dtype=numpy.float64)
747 elif isinstance(arg1,Symbol):
748 pass
749 else:
750 raise TypeError("function: Unknown type of second argument.")
751 else:
752 raise TypeError("function: Unknown type of first argument.")
753
754 return arg0,arg1
755
756 def matchShape(arg0,arg1):
757 """
758 Returns a representation of ``arg0`` and ``arg1`` which have the same shape.
759
760 :param arg0: first argument
761 :type arg0: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, `Symbol`
762 :param arg1: second argument
763 :type arg1: ``numpy.ndarray``,`escript.Data`,``float``, ``int``, `Symbol`
764 :return: ``arg0`` and ``arg1`` where copies are returned when the shape has
765 to be changed
766 :rtype: ``tuple``
767 """
768 sh=commonShape(arg0,arg1)
769 sh0=getShape(arg0)
770 sh1=getShape(arg1)
771 if len(sh0)<len(sh):
772 return outer(arg0,numpy.ones(sh[len(sh0):],numpy.float64)),arg1
773 elif len(sh1)<len(sh):
774 return arg0,outer(arg1,numpy.ones(sh[len(sh1):],numpy.float64))
775 else:
776 return arg0,arg1
777
778 def log10(arg):
779 """
780 Returns base-10 logarithm of argument ``arg``.
781
782 :param arg: argument
783 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
784 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
785 on the type of ``arg``
786 :raise TypeError: if the type of the argument is not expected
787 """
788 if isinstance(arg,numpy.ndarray):
789 return numpy.log10(arg)
790 elif isinstance(arg,escript.Data):
791 return arg._log10()
792 elif isinstance(arg,float):
793 return math.log10(arg)
794 elif isinstance(arg,int):
795 return math.log10(float(arg))
796 elif isinstance(arg,Symbol):
797 return arg.applyfunc(symfn.log10)
798 else:
799 raise TypeError("log10: Unknown argument type.")
800
801 def wherePositive(arg):
802 """
803 Returns mask of positive values of argument ``arg``.
804
805 :param arg: argument
806 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
807 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
808 on the type of ``arg``
809 :raise TypeError: if the type of the argument is not expected
810 """
811 if isinstance(arg,numpy.ndarray):
812 out=numpy.greater(arg,numpy.zeros(arg.shape,numpy.float64))*1.
813 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
814 return out
815 elif isinstance(arg,escript.Data):
816 return arg._wherePositive()
817 elif isinstance(arg,float):
818 if arg>0:
819 return 1.
820 else:
821 return 0.
822 elif isinstance(arg,int):
823 if arg>0:
824 return 1.
825 else:
826 return 0.
827 elif isinstance(arg,Symbol):
828 return arg.applyfunc(symfn.wherePositive)
829 else:
830 raise TypeError("wherePositive: Unknown argument type.")
831
832 def whereNegative(arg):
833 """
834 Returns mask of negative values of argument ``arg``.
835
836 :param arg: argument
837 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
838 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
839 on the type of ``arg``
840 :raise TypeError: if the type of the argument is not expected
841 """
842 if isinstance(arg,numpy.ndarray):
843 out=numpy.less(arg,numpy.zeros(arg.shape,numpy.float64))*1.
844 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
845 return out
846 elif isinstance(arg,escript.Data):
847 return arg._whereNegative()
848 elif isinstance(arg,float):
849 if arg<0:
850 return 1.
851 else:
852 return 0.
853 elif isinstance(arg,int):
854 if arg<0:
855 return 1.
856 else:
857 return 0.
858 elif isinstance(arg,Symbol):
859 return arg.applyfunc(symfn.whereNegative)
860 else:
861 raise TypeError("whereNegative: Unknown argument type.")
862
863 def whereNonNegative(arg):
864 """
865 Returns mask of non-negative values of argument ``arg``.
866
867 :param arg: argument
868 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
869 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
870 on the type of ``arg``
871 :raise TypeError: if the type of the argument is not expected
872 """
873 if isinstance(arg,numpy.ndarray):
874 out=numpy.greater_equal(arg,numpy.zeros(arg.shape,numpy.float64))*1.
875 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
876 return out
877 elif isinstance(arg,escript.Data):
878 return arg._whereNonNegative()
879 elif isinstance(arg,float):
880 if arg<0:
881 return 0.
882 else:
883 return 1.
884 elif isinstance(arg,int):
885 if arg<0:
886 return 0.
887 else:
888 return 1.
889 elif isinstance(arg,Symbol):
890 return arg.applyfunc(symfn.whereNonNegative)
891 else:
892 raise TypeError("whereNonNegative: Unknown argument type.")
893
894 def whereNonPositive(arg):
895 """
896 Returns mask of non-positive values of argument ``arg``.
897
898 :param arg: argument
899 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
900 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
901 on the type of ``arg``
902 :raise TypeError: if the type of the argument is not expected
903 """
904 if isinstance(arg,numpy.ndarray):
905 out=numpy.less_equal(arg,numpy.zeros(arg.shape,numpy.float64))*1.
906 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
907 return out
908 elif isinstance(arg,escript.Data):
909 return arg._whereNonPositive()
910 elif isinstance(arg,float):
911 if arg>0:
912 return 0.
913 else:
914 return 1.
915 elif isinstance(arg,int):
916 if arg>0:
917 return 0.
918 else:
919 return 1.
920 elif isinstance(arg,Symbol):
921 return arg.applyfunc(symfn.whereNonPositive)
922 else:
923 raise TypeError("whereNonPositive: Unknown argument type.")
924
925 def whereZero(arg,tol=None,rtol=math.sqrt(EPSILON)):
926 """
927 Returns mask of zero entries of argument ``arg``.
928
929 :param arg: argument
930 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
931 :param tol: absolute tolerance. Values with absolute value less than tol are accepted
932 as zero. If ``tol`` is not present ``rtol``*```Lsup` (arg)`` is used.
933 :type tol: ``float``
934 :param rtol: relative tolerance used to define the absolute tolerance if ``tol`` is not present.
935 :type rtol: non-negative ``float``
936 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
937 on the type of ``arg``
938 :raise ValueError: if ``rtol`` is non-negative.
939 :raise TypeError: if the type of the argument is not expected
940 """
941 if tol == None and not isinstance(arg, Symbol):
942 if rtol<0: raise ValueError("rtol must be non-negative.")
943 tol = Lsup(arg)*rtol
944 if isinstance(arg,numpy.ndarray):
945 out=numpy.less_equal(abs(arg)-tol,numpy.zeros(arg.shape,numpy.float64))*1.
946 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
947 return out
948 elif isinstance(arg,escript.Data):
949 return arg._whereZero(tol)
950 elif isinstance(arg,float):
951 if abs(arg)<=tol:
952 return 1.
953 else:
954 return 0.
955 elif isinstance(arg,int):
956 if abs(float(arg))<=tol:
957 return 1.
958 else:
959 return 0.
960 elif isinstance(arg,Symbol):
961 return arg.applyfunc(symfn.whereZero)
962 else:
963 raise TypeError("whereZero: Unknown argument type.")
964
965 def whereNonZero(arg,tol=0.):
966 """
967 Returns mask of values different from zero of argument ``arg``.
968
969 :param arg: argument
970 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
971 :param tol: absolute tolerance. Values with absolute value less than tol are accepted
972 as zero. If ``tol`` is not present ``rtol``*```Lsup` (arg)`` is used.
973 :type tol: ``float``
974 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
975 on the type of ``arg``
976 :raise ValueError: if ``rtol`` is non-negative.
977 :raise TypeError: if the type of the argument is not expected
978 """
979 if tol == None:
980 if rtol<=0: raise ValueError("rtol must be non-negative.")
981 tol = Lsup(arg)*rtol
982 if isinstance(arg,numpy.ndarray):
983 out=numpy.greater(abs(arg)-tol,numpy.zeros(arg.shape,numpy.float64))*1.
984 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
985 return out
986 elif isinstance(arg,escript.Data):
987 return arg._whereNonZero(tol)
988 elif isinstance(arg,float):
989 if abs(arg)>tol:
990 return 1.
991 else:
992 return 0.
993 elif isinstance(arg,int):
994 if abs(float(arg))>tol:
995 return 1.
996 else:
997 return 0.
998 elif isinstance(arg,Symbol):
999 return arg.applyfunc(symfn.whereNonZero)
1000 else:
1001 raise TypeError("whereNonZero: Unknown argument type.")
1002
1003 def Abs(arg):
1004 """
1005 Returns the absolute value of argument ``arg``.
1006
1007 :param arg: argument
1008 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1009 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1010 on the type of ``arg``
1011 :raise TypeError: if the type of the argument is not expected
1012 """
1013 if isinstance(arg,Symbol):
1014 return arg.applyfunc(symfn.abs)
1015 else:
1016 return abs(arg)
1017
1018 def erf(arg):
1019 """
1020 Returns the error function *erf* of argument ``arg``.
1021
1022 :param arg: argument
1023 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1024 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1025 on the type of ``arg``
1026 :raise TypeError: if the type of the argument is not expected
1027 """
1028 if isinstance(arg,escript.Data):
1029 return arg._erf()
1030 elif isinstance(arg,Symbol):
1031 return arg.applyfunc(symfn.erf)
1032 else:
1033 raise TypeError("erf: Unknown argument type.")
1034
1035 def sin(arg):
1036 """
1037 Returns sine of argument ``arg``.
1038
1039 :param arg: argument
1040 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1041 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1042 on the type of ``arg``
1043 :raise TypeError: if the type of the argument is not expected
1044 """
1045 if isinstance(arg,numpy.ndarray):
1046 return numpy.sin(arg)
1047 elif isinstance(arg,escript.Data):
1048 return arg._sin()
1049 elif isinstance(arg,float):
1050 return math.sin(arg)
1051 elif isinstance(arg,int):
1052 return math.sin(arg)
1053 elif isinstance(arg,Symbol):
1054 return arg.applyfunc(symfn.sin)
1055 else:
1056 raise TypeError("sin: Unknown argument type.")
1057
1058 def cos(arg):
1059 """
1060 Returns cosine of argument ``arg``.
1061
1062 :param arg: argument
1063 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1064 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1065 on the type of ``arg``
1066 :raise TypeError: if the type of the argument is not expected
1067 """
1068 if isinstance(arg,numpy.ndarray):
1069 return numpy.cos(arg)
1070 elif isinstance(arg,escript.Data):
1071 return arg._cos()
1072 elif isinstance(arg,float):
1073 return math.cos(arg)
1074 elif isinstance(arg,int):
1075 return math.cos(arg)
1076 elif isinstance(arg,Symbol):
1077 return arg.applyfunc(symfn.cos)
1078 else:
1079 raise TypeError("cos: Unknown argument type.")
1080
1081 def tan(arg):
1082 """
1083 Returns tangent of argument ``arg``.
1084
1085 :param arg: argument
1086 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1087 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1088 on the type of ``arg``
1089 :raise TypeError: if the type of the argument is not expected
1090 """
1091 if isinstance(arg,numpy.ndarray):
1092 return numpy.tan(arg)
1093 elif isinstance(arg,escript.Data):
1094 return arg._tan()
1095 elif isinstance(arg,float):
1096 return math.tan(arg)
1097 elif isinstance(arg,int):
1098 return math.tan(arg)
1099 elif isinstance(arg,Symbol):
1100 return arg.applyfunc(symfn.tan)
1101 else:
1102 raise TypeError("tan: Unknown argument type.")
1103
1104 def asin(arg):
1105 """
1106 Returns the inverse sine of argument ``arg``.
1107
1108 :param arg: argument
1109 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1110 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1111 on the type of ``arg``
1112 :raise TypeError: if the type of the argument is not expected
1113 """
1114 if isinstance(arg,numpy.ndarray):
1115 return numpy.arcsin(arg)
1116 elif isinstance(arg,escript.Data):
1117 return arg._asin()
1118 elif isinstance(arg,float):
1119 return math.asin(arg)
1120 elif isinstance(arg,int):
1121 return math.asin(arg)
1122 elif isinstance(arg,Symbol):
1123 return arg.applyfunc(symfn.asin)
1124 else:
1125 raise TypeError("asin: Unknown argument type.")
1126
1127 def acos(arg):
1128 """
1129 Returns the inverse cosine of argument ``arg``.
1130
1131 :param arg: argument
1132 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1133 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1134 on the type of ``arg``
1135 :raise TypeError: if the type of the argument is not expected
1136 """
1137 if isinstance(arg,numpy.ndarray):
1138 return numpy.arccos(arg)
1139 elif isinstance(arg,escript.Data):
1140 return arg._acos()
1141 elif isinstance(arg,float):
1142 return math.acos(arg)
1143 elif isinstance(arg,int):
1144 return math.acos(arg)
1145 elif isinstance(arg,Symbol):
1146 return arg.applyfunc(symfn.acos)
1147 else:
1148 raise TypeError("acos: Unknown argument type.")
1149
1150 def atan(arg):
1151 """
1152 Returns inverse tangent of argument ``arg``.
1153
1154 :param arg: argument
1155 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1156 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1157 on the type of ``arg``
1158 :raise TypeError: if the type of the argument is not expected
1159 """
1160 if isinstance(arg,numpy.ndarray):
1161 return numpy.arctan(arg)
1162 elif isinstance(arg,escript.Data):
1163 return arg._atan()
1164 elif isinstance(arg,float):
1165 return math.atan(arg)
1166 elif isinstance(arg,int):
1167 return math.atan(arg)
1168 elif isinstance(arg,Symbol):
1169 return arg.applyfunc(symfn.atan)
1170 else:
1171 raise TypeError("atan: Unknown argument type.")
1172
1173 def sinh(arg):
1174 """
1175 Returns the hyperbolic sine of argument ``arg``.
1176
1177 :param arg: argument
1178 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1179 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1180 on the type of ``arg``
1181 :raise TypeError: if the type of the argument is not expected
1182 """
1183 if isinstance(arg,numpy.ndarray):
1184 return numpy.sinh(arg)
1185 elif isinstance(arg,escript.Data):
1186 return arg._sinh()
1187 elif isinstance(arg,float):
1188 return math.sinh(arg)
1189 elif isinstance(arg,int):
1190 return math.sinh(arg)
1191 elif isinstance(arg,Symbol):
1192 return arg.applyfunc(symfn.sinh)
1193 else:
1194 raise TypeError("sinh: Unknown argument type.")
1195
1196 def cosh(arg):
1197 """
1198 Returns the hyperbolic cosine of argument ``arg``.
1199
1200 :param arg: argument
1201 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1202 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1203 on the type of ``arg``
1204 :raise TypeError: if the type of the argument is not expected
1205 """
1206 if isinstance(arg,numpy.ndarray):
1207 return numpy.cosh(arg)
1208 elif isinstance(arg,escript.Data):
1209 return arg._cosh()
1210 elif isinstance(arg,float):
1211 return math.cosh(arg)
1212 elif isinstance(arg,int):
1213 return math.cosh(arg)
1214 elif isinstance(arg,Symbol):
1215 return arg.applyfunc(symfn.cosh)
1216 else:
1217 raise TypeError("cosh: Unknown argument type.")
1218
1219 def tanh(arg):
1220 """
1221 Returns the hyperbolic tangent of argument ``arg``.
1222
1223 :param arg: argument
1224 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1225 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1226 on the type of ``arg``
1227 :raise TypeError: if the type of the argument is not expected
1228 """
1229 if isinstance(arg,numpy.ndarray):
1230 return numpy.tanh(arg)
1231 elif isinstance(arg,escript.Data):
1232 return arg._tanh()
1233 elif isinstance(arg,float):
1234 return math.tanh(arg)
1235 elif isinstance(arg,int):
1236 return math.tanh(arg)
1237 elif isinstance(arg,Symbol):
1238 return arg.applyfunc(symfn.tanh)
1239 else:
1240 raise TypeError("tanh: Unknown argument type.")
1241
1242 def asinh(arg):
1243 """
1244 Returns the inverse hyperbolic sine of argument ``arg``.
1245
1246 :param arg: argument
1247 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1248 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1249 on the type of ``arg``
1250 :raise TypeError: if the type of the argument is not expected
1251 """
1252 if isinstance(arg,numpy.ndarray):
1253 return numpy.arcsinh(arg)
1254 elif isinstance(arg,escript.Data):
1255 return arg._asinh()
1256 elif isinstance(arg,float):
1257 return numpy.arcsinh(arg)
1258 elif isinstance(arg,int):
1259 return numpy.arcsinh(float(arg))
1260 elif isinstance(arg,Symbol):
1261 return arg.applyfunc(symfn.asinh)
1262 else:
1263 raise TypeError("asinh: Unknown argument type.")
1264
1265 def acosh(arg):
1266 """
1267 Returns the inverse hyperbolic cosine of argument ``arg``.
1268
1269 :param arg: argument
1270 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1271 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1272 on the type of ``arg``
1273 :raise TypeError: if the type of the argument is not expected
1274 """
1275 if isinstance(arg,numpy.ndarray):
1276 return numpy.arccosh(arg)
1277 elif isinstance(arg,escript.Data):
1278 return arg._acosh()
1279 elif isinstance(arg,float):
1280 return numpy.arccosh(arg)
1281 elif isinstance(arg,int):
1282 return numpy.arccosh(float(arg))
1283 elif isinstance(arg,Symbol):
1284 return arg.applyfunc(symfn.acosh)
1285 else:
1286 raise TypeError("acosh: Unknown argument type.")
1287
1288 def atanh(arg):
1289 """
1290 Returns the inverse hyperbolic tangent of argument ``arg``.
1291
1292 :param arg: argument
1293 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1294 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1295 on the type of ``arg``
1296 :raise TypeError: if the type of the argument is not expected
1297 """
1298 if isinstance(arg,numpy.ndarray):
1299 return numpy.arctanh(arg)
1300 elif isinstance(arg,escript.Data):
1301 return arg._atanh()
1302 elif isinstance(arg,float):
1303 return numpy.arctanh(arg)
1304 elif isinstance(arg,int):
1305 return numpy.arctanh(float(arg))
1306 elif isinstance(arg,Symbol):
1307 return arg.applyfunc(symfn.atanh)
1308 else:
1309 raise TypeError("atanh: Unknown argument type.")
1310
1311 def exp(arg):
1312 """
1313 Returns *e* to the power of argument ``arg``.
1314
1315 :param arg: argument
1316 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1317 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1318 on the type of arg
1319 :raise TypeError: if the type of the argument is not expected
1320 """
1321 if isinstance(arg,numpy.ndarray):
1322 return numpy.exp(arg)
1323 elif isinstance(arg,escript.Data):
1324 return arg._exp()
1325 elif isinstance(arg,float):
1326 return math.exp(arg)
1327 elif isinstance(arg,int):
1328 return math.exp(arg)
1329 elif isinstance(arg,Symbol):
1330 return arg.applyfunc(symfn.exp)
1331 else:
1332 raise TypeError("exp: Unknown argument type.")
1333
1334 def sqrt(arg):
1335 """
1336 Returns the square root of argument ``arg``.
1337
1338 :param arg: argument
1339 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1340 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1341 depending on the type of ``arg``
1342 :raise TypeError: if the type of the argument is not expected
1343 """
1344 if isinstance(arg,numpy.ndarray):
1345 return numpy.sqrt(arg)
1346 elif isinstance(arg,escript.Data):
1347 return arg._sqrt()
1348 elif isinstance(arg,float):
1349 return math.sqrt(arg)
1350 elif isinstance(arg,int):
1351 return math.sqrt(arg)
1352 elif isinstance(arg,Symbol):
1353 return arg.applyfunc(symfn.sqrt)
1354 else:
1355 raise TypeError("sqrt: Unknown argument type.")
1356
1357 def log(arg):
1358 """
1359 Returns the natural logarithm of argument ``arg``.
1360
1361 :param arg: argument
1362 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``.
1363 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1364 on the type of ``arg``
1365 :raise TypeError: if the type of the argument is not expected
1366 """
1367 if isinstance(arg,numpy.ndarray):
1368 return numpy.log(arg)
1369 elif isinstance(arg,escript.Data):
1370 return arg._log()
1371 elif isinstance(arg,float):
1372 return math.log(arg)
1373 elif isinstance(arg,int):
1374 return math.log(arg)
1375 elif isinstance(arg,Symbol):
1376 return arg.applyfunc(symfn.log)
1377 else:
1378 raise TypeError("log: Unknown argument type.")
1379
1380 def sign(arg):
1381 """
1382 Returns the sign of argument ``arg``.
1383
1384 :param arg: argument
1385 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1386 :rtype: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray`` depending
1387 on the type of ``arg``
1388 :raise TypeError: if the type of the argument is not expected
1389 """
1390 if isinstance(arg,numpy.ndarray):
1391 return wherePositive(arg)-whereNegative(arg)
1392 elif isinstance(arg,escript.Data):
1393 return arg._sign()
1394 elif isinstance(arg,float):
1395 if arg>0:
1396 return 1.
1397 elif arg<0:
1398 return -1.
1399 else:
1400 return 0.
1401 elif isinstance(arg,int):
1402 if float(arg)>0:
1403 return 1.
1404 elif float(arg)<0:
1405 return -1.
1406 else:
1407 return 0.
1408 elif isinstance(arg,Symbol):
1409 return arg.applyfunc(symfn.sign)
1410 else:
1411 raise TypeError("sign: Unknown argument type.")
1412
1413 def minval(arg):
1414 """
1415 Returns the minimum value over all components of ``arg`` at each data point.
1416
1417 :param arg: argument
1418 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1419 :rtype: ``float``, `escript.Data`, `Symbol` depending on the type of ``arg``
1420 :raise TypeError: if the type of the argument is not expected
1421 """
1422 if isinstance(arg,numpy.ndarray):
1423 if arg.ndim==0:
1424 return float(arg)
1425 else:
1426 return arg.min()
1427 elif isinstance(arg,escript.Data):
1428 return arg._minval()
1429 elif isinstance(arg,float):
1430 return arg
1431 elif isinstance(arg,int):
1432 return float(arg)
1433 elif isinstance(arg,Symbol):
1434 return symfn.minval(arg)
1435 else:
1436 raise TypeError("minval: Unknown argument type.")
1437
1438 def maxval(arg):
1439 """
1440 Returns the maximum value over all components of ``arg`` at each data point.
1441
1442 :param arg: argument
1443 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1444 :rtype: ``float``, `escript.Data`, `Symbol` depending on the type of ``arg``
1445 :raise TypeError: if the type of the argument is not expected
1446 """
1447 if isinstance(arg,numpy.ndarray):
1448 if arg.ndim==0:
1449 return float(arg)
1450 else:
1451 return arg.max()
1452 elif isinstance(arg,escript.Data):
1453 return arg._maxval()
1454 elif isinstance(arg,float):
1455 return arg
1456 elif isinstance(arg,int):
1457 return float(arg)
1458 elif isinstance(arg,Symbol):
1459 return symfn.maxval(arg)
1460 else:
1461 raise TypeError("maxval: Unknown argument type.")
1462
1463 def length(arg):
1464 """
1465 Returns the length (Euclidean norm) of argument ``arg`` at each data point.
1466
1467 :param arg: argument
1468 :type arg: ``float``, `escript.Data`, `Symbol`, ``numpy.ndarray``
1469 :rtype: ``float``, `escript.Data`, `Symbol` depending on the type of ``arg``
1470 """
1471 return sqrt(inner(arg,arg))
1472
1473 def trace(arg,axis_offset=0):
1474 """
1475 Returns the trace of ``arg`` which is the sum of ``arg[k,k]`` over k.
1476
1477 :param arg: argument
1478 :type arg: `escript.Data`, `Symbol`, ``numpy.ndarray``
1479 :param axis_offset: ``axis_offset`` to components to sum over. ``axis_offset``
1480 must be non-negative and less than the rank of ``arg`` +1.
1481 The dimensions of component ``axis_offset`` and
1482 axis_offset+1 must be equal.
1483 :type axis_offset: ``int``
1484 :return: trace of arg. The rank of the returned object is rank of ``arg``
1485 minus 2.
1486 :rtype: `escript.Data`, `Symbol` or ``numpy.ndarray`` depending on the type
1487 of ``arg``
1488 """
1489 if isinstance(arg,numpy.ndarray):
1490 sh=arg.shape
1491 if len(sh)<2:
1492 raise ValueError("rank of argument must be greater than 1")
1493 if axis_offset<0 or axis_offset>len(sh)-2:
1494 raise ValueError("axis_offset must be between 0 and %d"%(len(sh)-2))
1495 s1=1
1496 for i in range(axis_offset): s1*=sh[i]
1497 s2=1
1498 for i in range(axis_offset+2,len(sh)): s2*=sh[i]
1499 if not sh[axis_offset] == sh[axis_offset+1]:
1500 raise ValueError("dimensions of component %d and %d must match."%(axis_offset,axis_offset+1))
1501 arg_reshaped=numpy.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
1502 out=numpy.zeros([s1,s2],numpy.float64)
1503 for i1 in range(s1):
1504 for i2 in range(s2):
1505 for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
1506 out.resize(sh[:axis_offset]+sh[axis_offset+2:])
1507 return out
1508 elif isinstance(arg,escript.Data):
1509 if arg.getRank()<2:
1510 raise ValueError("rank of argument must be greater than 1")
1511 if axis_offset<0 or axis_offset>arg.getRank()-2:
1512 raise ValueError("axis_offset must be between 0 and %d"%(arg.getRank()-2))
1513 s=list(arg.getShape())
1514 if not s[axis_offset] == s[axis_offset+1]:
1515 raise ValueError("dimensions of component %d and %d must match."%(axis_offset,axis_offset+1))
1516 return arg._trace(axis_offset)
1517 elif isinstance(arg,Symbol):
1518 if arg.getRank()<2:
1519 raise ValueError("rank of argument must be greater than 1")
1520 if axis_offset<0 or axis_offset>arg.getRank()-2:
1521 raise ValueError("axis_offset must be between 0 and %d"%(arg.getRank()-2))
1522 s=list(arg.getShape())
1523 if not s[axis_offset] == s[axis_offset+1]:
1524 raise ValueError("dimensions of component %d and %d must match."%(axis_offset,axis_offset+1))
1525 return arg.trace(axis_offset)
1526 elif isinstance(arg,float):
1527 raise TypeError,"illegal argument type float."
1528 elif isinstance(arg,int):
1529 raise TypeError,"illegal argument type int."
1530 else:
1531 raise TypeError,"Unknown argument type."
1532
1533 def transpose(arg,axis_offset=None):
1534 """
1535 Returns the transpose of ``arg`` by swapping the first ``axis_offset`` and
1536 the last ``rank-axis_offset`` components.
1537
1538 :param arg: argument
1539 :type arg: `escript.Data`, `Symbol`, ``numpy.ndarray``, ``float``, ``int``
1540 :param axis_offset: the first ``axis_offset`` components are swapped with the
1541 rest. ``axis_offset`` must be non-negative and less or
1542 equal to the rank of ``arg``. If ``axis_offset`` is not
1543 present ``int(r/2)`` where r is the rank of ``arg`` is
1544 used.
1545 :type axis_offset: ``int``
1546 :return: transpose of ``arg``
1547 :rtype: `escript.Data`, `Symbol`, ``numpy.ndarray``, ``float``, ``int``
1548 depending on the type of ``arg``
1549 """
1550 if isinstance(arg,numpy.ndarray):
1551 if axis_offset==None: axis_offset=int(arg.ndim/2)
1552 return numpy.transpose(arg,axes=list(range(axis_offset,arg.ndim))+list(range(0,axis_offset)))
1553 elif isinstance(arg,escript.Data):
1554 r=arg.getRank()
1555 if axis_offset==None: axis_offset=int(r/2)
1556 if axis_offset<0 or axis_offset>r:
1557 raise ValueError("axis_offset must be between 0 and %s"%r)
1558 return arg._transpose(axis_offset)
1559 elif isinstance(arg,float):
1560 if not ( axis_offset==0 or axis_offset==None):
1561 raise ValueError("axis_offset must be 0 for float argument")
1562 return arg
1563 elif isinstance(arg,int):
1564 if not ( axis_offset==0 or axis_offset==None):
1565 raise ValueError("axis_offset must be 0 for int argument")
1566 return float(arg)
1567 elif isinstance(arg,Symbol):
1568 r=arg.getRank()
1569 if axis_offset==None: axis_offset=int(r/2)
1570 if axis_offset<0 or axis_offset>r:
1571 raise ValueError,"axis_offset must be between 0 and %s"%r
1572 return arg.transpose(axis_offset)
1573 else:
1574 raise TypeError("Unknown argument type.")
1575
1576 def swap_axes(arg,axis0=0,axis1=1):
1577 """
1578 Returns the swap of ``arg`` by swapping the components ``axis0`` and ``axis1``.
1579
1580 :param arg: argument
1581 :type arg: `escript.Data`, `Symbol`, ``numpy.ndarray``
1582 :param axis0: first axis. ``axis0`` must be non-negative and less than the
1583 rank of ``arg``.
1584 :type axis0: ``int``
1585 :param axis1: second axis. ``axis1`` must be non-negative and less than the
1586 rank of ``arg``.
1587 :type axis1: ``int``
1588 :return: ``arg`` with swapped components
1589 :rtype: `escript.Data`, `Symbol` or ``numpy.ndarray`` depending on the type
1590 of ``arg``
1591 """
1592 if axis0 > axis1:
1593 axis0,axis1=axis1,axis0
1594 if isinstance(arg,numpy.ndarray):
1595 return numpy.swapaxes(arg,axis0,axis1)
1596 elif isinstance(arg,escript.Data):
1597 return arg._swap_axes(axis0,axis1)
1598 elif isinstance(arg,Symbol):
1599 return arg.swap_axes(axis0,axis1)
1600 elif isinstance(arg,float):
1601 raise TypeError("float argument is not supported.")
1602 elif isinstance(arg,int):
1603 raise TypeError("int argument is not supported.")
1604 else:
1605 raise TypeError("Unknown argument type.")
1606
1607 def symmetric(arg):
1608 """
1609 Returns the symmetric part of the square matrix ``arg``. That is,
1610 *(arg+transpose(arg))/2*.
1611
1612 :param arg: input matrix. Must have rank 2 or 4 and be square.
1613 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1614 :return: symmetric part of ``arg``
1615 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1616 """
1617 if isinstance(arg,numpy.ndarray):
1618 if arg.ndim==2:
1619 if not (arg.shape[0]==arg.shape[1]):
1620 raise ValueError("argument must be square.")
1621 elif arg.ndim==4:
1622 if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
1623 raise ValueError("argument must be square.")
1624 else:
1625 raise ValueError("rank 2 or 4 is required.")
1626 return (arg+transpose(arg))/2
1627 elif isinstance(arg,escript.Data):
1628 if arg.getRank()==2:
1629 if not (arg.getShape()[0]==arg.getShape()[1]):
1630 raise ValueError("argument must be square.")
1631 return arg._symmetric()
1632 elif arg.getRank()==4:
1633 if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
1634 raise ValueError("argument must be square.")
1635 return arg._symmetric()
1636 else:
1637 raise ValueError("rank 2 or 4 is required.")
1638 elif isinstance(arg, Symbol):
1639 if arg.getRank()==2:
1640 if arg.getShape()[0]!=arg.getShape()[1]:
1641 raise ValueError("symmetric: argument must be square.")
1642 elif arg.getRank()==4:
1643 if arg.getShape()[0]!=arg.getShape()[2] or arg.getShape()[1]!=arg.getShape()[3]:
1644 raise ValueError("symmetric: argument must be square.")
1645 else:
1646 raise ValueError("symmetric: rank 2 or 4 is required.")
1647 return (arg+transpose(arg))/2
1648 elif isinstance(arg,float):
1649 return arg
1650 elif isinstance(arg,int):
1651 return float(arg)
1652 else:
1653 raise TypeError("symmetric: Unknown argument type.")
1654
1655 def nonsymmetric(arg):
1656 """
1657 Returns the non-symmetric part of the square matrix ``arg``. That is,
1658 *(arg-transpose(arg))/2*.
1659
1660 :param arg: input matrix. Must have rank 2 or 4 and be square.
1661 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1662 :return: non-symmetric part of ``arg``
1663 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1664 """
1665 if isinstance(arg,numpy.ndarray):
1666 if arg.ndim==2:
1667 if not (arg.shape[0]==arg.shape[1]):
1668 raise ValueError("nonsymmetric: argument must be square.")
1669 elif arg.ndim==4:
1670 if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
1671 raise ValueError("nonsymmetric: argument must be square.")
1672 else:
1673 raise ValueError("nonsymmetric: rank 2 or 4 is required.")
1674 return (arg-transpose(arg))/2
1675 elif isinstance(arg,escript.Data):
1676 if arg.getRank()==2:
1677 if not (arg.getShape()[0]==arg.getShape()[1]):
1678 raise ValueError("argument must be square.")
1679 return arg._nonsymmetric()
1680 elif arg.getRank()==4:
1681 if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
1682 raise ValueError("argument must be square.")
1683 return arg._nonsymmetric()
1684 else:
1685 raise ValueError("rank 2 or 4 is required.")
1686 elif isinstance(arg, Symbol):
1687 if arg.getRank()==2:
1688 if arg.getShape()[0]!=arg.getShape()[1]:
1689 raise ValueError("nonsymmetric: argument must be square.")
1690 elif arg.getRank()==4:
1691 if arg.getShape()[0]!=arg.getShape()[2] or arg.getShape()[1]!=arg.getShape()[3]:
1692 raise ValueError("nonsymmetric: argument must be square.")
1693 else:
1694 raise ValueError("nonsymmetric: rank 2 or 4 is required.")
1695 return (arg-transpose(arg))/2
1696 elif isinstance(arg,float):
1697 return arg
1698 elif isinstance(arg,int):
1699 return float(arg)
1700 else:
1701 raise TypeError("nonsymmetric: Unknown argument type.")
1702
1703 def inverse(arg):
1704 """
1705 Returns the inverse of the square matrix ``arg``.
1706
1707 :param arg: square matrix. Must have rank 2 and the first and second
1708 dimension must be equal.
1709 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1710 :return: inverse of the argument. ``matrix_mult(inverse(arg),arg)`` will be
1711 almost equal to ``kronecker(arg.getShape()[0])``
1712 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1713 :note: for `escript.Data` objects the dimension is restricted to 3.
1714 """
1715 import numpy.linalg
1716 if isinstance(arg,numpy.ndarray):
1717 return numpy.linalg.tensorinv(arg,ind=1)
1718 elif isinstance(arg,escript.Data):
1719 return escript_inverse(arg)
1720 elif isinstance(arg,float):
1721 return 1./arg
1722 elif isinstance(arg,int):
1723 return 1./float(arg)
1724 elif isinstance(arg,Symbol):
1725 return arg.inverse()
1726 else:
1727 raise TypeError("inverse: Unknown argument type.")
1728
1729 def escript_inverse(arg): # this should be escript._inverse and use LAPACK
1730 "arg is a Data object!"
1731 return arg._inverse()
1732
1733 #if not arg.getRank()==2:
1734 #raise ValueError,"escript_inverse: argument must have rank 2"
1735 #s=arg.getShape()
1736 #if not s[0] == s[1]:
1737 #raise ValueError,"escript_inverse: argument must be a square matrix."
1738 #out=escript.Data(0.,s,arg.getFunctionSpace())
1739 #if s[0]==1:
1740 #if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
1741 #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1742 #out[0,0]=1./arg[0,0]
1743 #elif s[0]==2:
1744 #A11=arg[0,0]
1745 #A12=arg[0,1]
1746 #A21=arg[1,0]
1747 #A22=arg[1,1]
1748 #D = A11*A22-A12*A21
1749 #if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
1750 #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1751 #D=1./D
1752 #out[0,0]= A22*D
1753 #out[1,0]=-A21*D
1754 #out[0,1]=-A12*D
1755 #out[1,1]= A11*D
1756 #elif s[0]==3:
1757 #A11=arg[0,0]
1758 #A21=arg[1,0]
1759 #A31=arg[2,0]
1760 #A12=arg[0,1]
1761 #A22=arg[1,1]
1762 #A32=arg[2,1]
1763 #A13=arg[0,2]
1764 #A23=arg[1,2]
1765 #A33=arg[2,2]
1766 #D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
1767 #if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
1768 #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1769 #D=1./D
1770 #out[0,0]=(A22*A33-A23*A32)*D
1771 #out[1,0]=(A31*A23-A21*A33)*D
1772 #out[2,0]=(A21*A32-A31*A22)*D
1773 #out[0,1]=(A13*A32-A12*A33)*D
1774 #out[1,1]=(A11*A33-A31*A13)*D
1775 #out[2,1]=(A12*A31-A11*A32)*D
1776 #out[0,2]=(A12*A23-A13*A22)*D
1777 #out[1,2]=(A13*A21-A11*A23)*D
1778 #out[2,2]=(A11*A22-A12*A21)*D
1779 #else:
1780 #raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
1781 #return out
1782
1783 def eigenvalues(arg):
1784 """
1785 Returns the eigenvalues of the square matrix ``arg``.
1786
1787 :param arg: square matrix. Must have rank 2 and the first and second
1788 dimension must be equal. It must also be symmetric, ie.
1789 ``transpose(arg)==arg`` (this is not checked).
1790 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`
1791 :return: the eigenvalues in increasing order
1792 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
1793 :note: for `escript.Data` and `Symbol` objects the dimension is
1794 restricted to 3.
1795 """
1796 if isinstance(arg,numpy.ndarray):
1797 out=numpy.linalg.eigvals((arg+numpy.transpose(arg))/2.)
1798 out.sort()
1799 return out
1800 elif isinstance(arg,escript.Data):
1801 return arg._eigenvalues()
1802 elif isinstance(arg,float):
1803 return arg
1804 elif isinstance(arg,int):
1805 return float(arg)
1806 elif isinstance(arg,Symbol):
1807 return symfn.eigenvalues(arg)
1808 else:
1809 raise TypeError("eigenvalues: Unknown argument type.")
1810
1811 def eigenvalues_and_eigenvectors(arg):
1812 """
1813 Returns the eigenvalues and eigenvectors of the square matrix ``arg``.
1814
1815 :param arg: square matrix. Must have rank 2 and the first and second
1816 dimension must be equal. It must also be symmetric, ie.
1817 ``transpose(arg)==arg`` (this is not checked).
1818 :type arg: `escript.Data`
1819 :return: the eigenvalues and eigenvectors. The eigenvalues are ordered by
1820 increasing value. The eigenvectors are orthogonal and normalized.
1821 If V are the eigenvectors then V[:,i] is the eigenvector
1822 corresponding to the i-th eigenvalue.
1823 :rtype: `tuple` of `escript.Data`
1824 :note: The dimension is restricted to 3.
1825 """
1826 if isinstance(arg,numpy.ndarray):
1827 raise TypeError("eigenvalues_and_eigenvectors does not support numpy.ndarray arguments")
1828 elif isinstance(arg,escript.Data):
1829 return arg._eigenvalues_and_eigenvectors()
1830 elif isinstance(arg,float):
1831 return (numpy.array([[arg]],numpy.float_),numpy.ones((1,1),numpy.float_))
1832 elif isinstance(arg,int):
1833 return (numpy.array([[arg]],numpy.float_),numpy.ones((1,1),numpy.float_))
1834 elif isinstance(arg,Symbol):
1835 return symfn.eigenvalues_and_eigenvectors(arg)
1836 else:
1837 raise TypeError("eigenvalues: Unknown argument type.")
1838
1839 def mult(arg0,arg1):
1840 """
1841 Product of ``arg0`` and ``arg1``.
1842
1843 :param arg0: first term
1844 :type arg0: `Symbol`, ``float``, ``int``, `escript.Data` or
1845 ``numpy.ndarray``
1846 :param arg1: second term
1847 :type arg1: `Symbol`, ``float``, ``int``, `escript.Data` or
1848 ``numpy.ndarray``
1849 :return: the product of ``arg0`` and ``arg1``
1850 :rtype: `Symbol`, ``float``, ``int``, `escript.Data` or
1851 ``numpy.ndarray``
1852 :note: The shape of both arguments is matched according to the rules
1853 used in `matchShape`.
1854 """
1855 args=matchShape(arg0,arg1)
1856 if testForZero(args[0]) or testForZero(args[1]):
1857 return numpy.zeros(getShape(args[0]),numpy.float64)
1858 else:
1859 if isinstance(args[0],numpy.ndarray):
1860 return args[1]*args[0]
1861 else:
1862 return args[0]*args[1]
1863
1864 def maximum(*args):
1865 """
1866 The maximum over arguments ``args``.
1867
1868 :param args: arguments
1869 :type args: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1870 ``float``
1871 :return: an object which in each entry gives the maximum of the
1872 corresponding values in ``args``
1873 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1874 ``float`` depending on the input
1875 """
1876 if max([isinstance(v,Symbol) for v in args]):
1877 return symfn.maximum(*args)
1878 out=None
1879 for a in args:
1880 if out==None:
1881 out=a*1.
1882 else:
1883 if isinstance(out,escript.Data) and isinstance(a,escript.Data):
1884 if out.getRank()==0 and a.getRank()>0:
1885 # We need to consider the case where we have scalars and
1886 # higher ranked objects mixed. If the scalar was first it will
1887 # get picked as the initial out and we have a problem, so we
1888 # swap the objects
1889 res=a.copy() #Deep copy of a
1890 res.copyWithMask(out,wherePositive(out-a))
1891 out=res
1892 else:
1893 out.copyWithMask(a,wherePositive(a-out))
1894 else:
1895 if isinstance(a, numpy.ndarray):
1896 diff=-out+a
1897 else:
1898 diff=a-out
1899 temp=mult(whereNonPositive(diff),out)+mult(wherePositive(diff),a)
1900 if isinstance(out,numpy.ndarray) and isinstance(a,numpy.ndarray):
1901 # we need to convert the result to an array
1902 temp=numpy.array(temp)
1903 out=temp
1904 return out
1905
1906 def minimum(*args):
1907 """
1908 The minimum over arguments ``args``.
1909
1910 :param args: arguments
1911 :type args: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1912 ``float``
1913 :return: an object which gives in each entry the minimum of the
1914 corresponding values in ``args``
1915 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1916 ``float`` depending on the input
1917 """
1918 if max([isinstance(v,Symbol) for v in args]):
1919 return symfn.minimum(*args)
1920 out=None
1921 for a in args:
1922 if out==None:
1923 out=a*1.
1924 else:
1925 if isinstance(out,escript.Data) and isinstance(a,escript.Data):
1926 if out.getRank()==0 and a.getRank()>0:
1927 # We need to consider the case where we have scalars and
1928 # higher ranked objects mixed. If the scalar was first it will
1929 # get picked as the initial out and we have a problem, so we
1930 # swap the objects
1931 res=a.copy() # Deep copy of a
1932 res.copyWithMask(out,whereNegative(out-a))
1933 out=res
1934 else:
1935 out.copyWithMask(a,whereNegative(a-out))
1936 else:
1937 if isinstance(a, numpy.ndarray):
1938 diff=-out+a
1939 else:
1940 diff=a-out
1941 #out=add(out,mult(whereNegative(diff),diff))
1942 temp=mult(whereNonNegative(diff),out)+mult(whereNegative(diff),a)
1943 if isinstance(out,numpy.ndarray) and isinstance(a,numpy.ndarray):
1944 # we need to convert the result to an array
1945 temp=numpy.array(temp)
1946 out=temp
1947 return out
1948
1949 def clip(arg,minval=None,maxval=None):
1950 """
1951 Cuts the values of ``arg`` between ``minval`` and ``maxval``.
1952
1953 :param arg: argument
1954 :type arg: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1955 ``float``
1956 :param minval: lower range. If None no lower range is applied
1957 :type minval: ``float`` or ``None``
1958 :param maxval: upper range. If None no upper range is applied
1959 :type maxval: ``float`` or ``None``
1960 :return: an object that contains all values from ``arg`` between ``minval``
1961 and ``maxval``
1962 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``int`` or
1963 ``float`` depending on the input
1964 :raise ValueError: if ``minval>maxval``
1965 """
1966 if isinstance(arg, Symbol):
1967 clip_item=lambda item: symfn.clip(item, minval, maxval)
1968 return arg.applyfunc(clip_item)
1969 if not minval==None and not maxval==None:
1970 if minval>maxval:
1971 raise ValueError("minval = %s must be less than maxval %s"%(minval,maxval))
1972 if minval == None:
1973 tmp=arg
1974 else:
1975 tmp=maximum(minval,arg)
1976 if maxval == None:
1977 return tmp
1978 else:
1979 return minimum(tmp,maxval)
1980
1981 def inner(arg0,arg1):
1982 """
1983 Inner product of the two arguments. The inner product is defined as:
1984
1985 C{out=Sigma_s arg0[s]*arg1[s]}
1986
1987 where s runs through ``arg0.Shape``.
1988
1989 ``arg0`` and ``arg1`` must have the same shape.
1990
1991 :param arg0: first argument
1992 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
1993 :param arg1: second argument
1994 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
1995 :return: the inner product of ``arg0`` and ``arg1`` at each data point
1996 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``
1997 depending on the input
1998 :raise ValueError: if the shapes of the arguments are not identical
1999 """
2000 sh0=getShape(arg0)
2001 sh1=getShape(arg1)
2002 if not sh0==sh1:
2003 raise ValueError("inner: shape of arguments does not match")
2004 return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
2005
2006 def outer(arg0,arg1):
2007 """
2008 The outer product of the two arguments. The outer product is defined as:
2009
2010 ``out[t,s]=arg0[t]*arg1[s]``
2011
2012 where
2013 - s runs through ``arg0.Shape``
2014 - t runs through ``arg1.Shape``
2015
2016 :param arg0: first argument
2017 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2018 :param arg1: second argument
2019 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2020 :return: the outer product of ``arg0`` and ``arg1`` at each data point
2021 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2022 """
2023 return generalTensorProduct(arg0,arg1,axis_offset=0)
2024
2025 def matrixmult(arg0,arg1):
2026 """
2027 See `matrix_mult`.
2028 """
2029 return matrix_mult(arg0,arg1)
2030
2031 def matrix_mult(arg0,arg1):
2032 """
2033 matrix-matrix or matrix-vector product of the two arguments.
2034
2035 C{out[s0]=Sigma_{r0} arg0[s0,r0]*arg1[r0]}
2036
2037 or
2038
2039 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[r0,s1]}
2040
2041 The second dimension of ``arg0`` and the first dimension of ``arg1`` must
2042 match.
2043
2044 :param arg0: first argument of rank 2
2045 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2046 :param arg1: second argument of at least rank 1
2047 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2048 :return: the matrix-matrix or matrix-vector product of ``arg0`` and ``arg1``
2049 at each data point
2050 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2051 :raise ValueError: if the shapes of the arguments are not appropriate
2052 """
2053 sh0=getShape(arg0)
2054 sh1=getShape(arg1)
2055 if not len(sh0)==2 :
2056 raise ValueError("first argument must have rank 2")
2057 if not len(sh1)==2 and not len(sh1)==1:
2058 raise ValueError("second argument must have rank 1 or 2")
2059 return generalTensorProduct(arg0,arg1,axis_offset=1)
2060
2061 def tensormult(arg0,arg1):
2062 """
2063 See `tensor_mult`.
2064 """
2065 return tensor_mult(arg0,arg1)
2066
2067 def tensor_mult(arg0,arg1):
2068 """
2069 The tensor product of the two arguments.
2070
2071 For ``arg0`` of rank 2 this is
2072
2073 C{out[s0]=Sigma_{r0} arg0[s0,r0]*arg1[r0]}
2074
2075 or
2076
2077 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[r0,s1]}
2078
2079 and for ``arg0`` of rank 4 this is
2080
2081 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]}
2082
2083 or
2084
2085 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]}
2086
2087 or
2088
2089 C{out[s0,s1]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]}
2090
2091 In the first case the second dimension of ``arg0`` and the last dimension of
2092 ``arg1`` must match and in the second case the two last dimensions of ``arg0``
2093 must match the two first dimensions of ``arg1``.
2094
2095 :param arg0: first argument of rank 2 or 4
2096 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2097 :param arg1: second argument of shape greater than 1 or 2 depending on the
2098 rank of ``arg0``
2099 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2100 :return: the tensor product of ``arg0`` and ``arg1`` at each data point
2101 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2102 """
2103 sh0=getShape(arg0)
2104 sh1=getShape(arg1)
2105 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2106 return generalTensorProduct(arg0,arg1,axis_offset=1)
2107 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2108 return generalTensorProduct(arg0,arg1,axis_offset=2)
2109 else:
2110 raise ValueError("tensor_mult: first argument must have rank 2 or 4")
2111
2112 def generalTensorProduct(arg0,arg1,axis_offset=0):
2113 """
2114 Generalized tensor product.
2115
2116 C{out[s,t]=Sigma_r arg0[s,r]*arg1[r,t]}
2117
2118 where
2119 - s runs through ``arg0.Shape[:arg0.ndim-axis_offset]``
2120 - r runs through ``arg1.Shape[:axis_offset]``
2121 - t runs through ``arg1.Shape[axis_offset:]``
2122
2123 :param arg0: first argument
2124 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2125 :param arg1: second argument
2126 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2127 :return: the general tensor product of ``arg0`` and ``arg1`` at each data
2128 point
2129 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2130 """
2131 if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2132 arg0,arg1=matchType(arg0,arg1)
2133 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data,
2134 # or one is a Symbol and the other either of the allowed types
2135 if isinstance(arg0,Symbol):
2136 sh0=arg0.getShape()
2137 sh1=getShape(arg1)
2138 if not sh0[arg0.getRank()-axis_offset:]==sh1[:axis_offset]:
2139 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))
2140 if isinstance(arg1,float):
2141 return arg0*arg1
2142 elif isinstance(arg1,numpy.ndarray) or isinstance(arg1, Symbol):
2143 return arg0.tensorProduct(arg1, axis_offset)
2144 elif isinstance(arg1, escript.Data):
2145 raise TypeError("tensor product of Symbol and Data not supported yet")
2146 elif isinstance(arg0,numpy.ndarray):
2147 if not arg0.shape[arg0.ndim-axis_offset:]==arg1.shape[:axis_offset]:
2148 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))
2149 arg0_c=arg0.copy()
2150 arg1_c=arg1.copy()
2151 sh0,sh1=arg0.shape,arg1.shape
2152 d0,d1,d01=1,1,1
2153 for i in sh0[:arg0.ndim-axis_offset]: d0*=i
2154 for i in sh1[axis_offset:]: d1*=i
2155 for i in sh1[:axis_offset]: d01*=i
2156 arg0_c.resize((d0,d01))
2157 arg1_c.resize((d01,d1))
2158 out=numpy.zeros((d0,d1),numpy.float64)
2159 for i0 in range(d0):
2160 for i1 in range(d1):
2161 out[i0,i1]=numpy.sum(arg0_c[i0,:]*arg1_c[:,i1])
2162 out.resize(sh0[:arg0.ndim-axis_offset]+sh1[axis_offset:])
2163 return out
2164 elif isinstance(arg0,escript.Data):
2165 if isinstance(arg1, Symbol):
2166 raise TypeError("tensor product of Data and Symbol not supported yet")
2167 return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2168 raise TypeError("generalTensorProduct: Unsupported argument type")
2169
2170 def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
2171 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2172 return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
2173
2174 def transposed_matrix_mult(arg0,arg1):
2175 """
2176 transposed(matrix)-matrix or transposed(matrix)-vector product of the two
2177 arguments.
2178
2179 C{out[s0]=Sigma_{r0} arg0[r0,s0]*arg1[r0]}
2180
2181 or
2182
2183 C{out[s0,s1]=Sigma_{r0} arg0[r0,s0]*arg1[r0,s1]}
2184
2185 The function call ``transposed_matrix_mult(arg0,arg1)`` is equivalent to
2186 ``matrix_mult(transpose(arg0),arg1)``.
2187
2188 The first dimension of ``arg0`` and ``arg1`` must match.
2189
2190 :param arg0: first argument of rank 2
2191 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2192 :param arg1: second argument of at least rank 1
2193 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2194 :return: the product of the transpose of ``arg0`` and ``arg1`` at each data
2195 point
2196 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2197 :raise ValueError: if the shapes of the arguments are not appropriate
2198 """
2199 sh0=getShape(arg0)
2200 sh1=getShape(arg1)
2201 if not len(sh0)==2 :
2202 raise ValueError("first argument must have rank 2")
2203 if not len(sh1)==2 and not len(sh1)==1:
2204 raise ValueError("second argument must have rank 1 or 2")
2205 return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
2206
2207 def transposed_tensor_mult(arg0,arg1):
2208 """
2209 The tensor product of the transpose of the first and the second argument.
2210
2211 For ``arg0`` of rank 2 this is
2212
2213 C{out[s0]=Sigma_{r0} arg0[r0,s0]*arg1[r0]}
2214
2215 or
2216
2217 C{out[s0,s1]=Sigma_{r0} arg0[r0,s0]*arg1[r0,s1]}
2218
2219 and for ``arg0`` of rank 4 this is
2220
2221 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]}
2222
2223 or
2224
2225 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]}
2226
2227 or
2228
2229 C{out[s0,s1]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]}
2230
2231 In the first case the first dimension of ``arg0`` and the first dimension of
2232 ``arg1`` must match and in the second case the two first dimensions of
2233 ``arg0`` must match the two first dimensions of ``arg1``.
2234
2235 The function call ``transposed_tensor_mult(arg0,arg1)`` is equivalent to
2236 ``tensor_mult(transpose(arg0),arg1)``.
2237
2238 :param arg0: first argument of rank 2 or 4
2239 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2240 :param arg1: second argument of shape greater of 1 or 2 depending on the
2241 rank of ``arg0``
2242 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2243 :return: the tensor product of transpose of arg0 and arg1 at each data point
2244 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2245 """
2246 sh0=getShape(arg0)
2247 sh1=getShape(arg1)
2248 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2249 return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
2250 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2251 return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
2252 else:
2253 raise ValueError("first argument must have rank 2 or 4")
2254
2255 def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
2256 """
2257 Generalized tensor product of transposed of ``arg0`` and ``arg1``.
2258
2259 C{out[s,t]=Sigma_r arg0[r,s]*arg1[r,t]}
2260
2261 where
2262 - s runs through ``arg0.Shape[axis_offset:]``
2263 - r runs through ``arg0.Shape[:axis_offset]``
2264 - t runs through ``arg1.Shape[axis_offset:]``
2265
2266 The function call ``generalTransposedTensorProduct(arg0,arg1,axis_offset)``
2267 is equivalent to
2268 ``generalTensorProduct(transpose(arg0,arg0.ndim-axis_offset),arg1,axis_offset)``.
2269
2270 :param arg0: first argument
2271 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2272 :param arg1: second argument
2273 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2274 :return: the general tensor product of ``transpose(arg0)`` and ``arg1`` at
2275 each data point
2276 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2277 """
2278 if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2279 arg0,arg1=matchType(arg0,arg1)
2280 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data,
2281 # or one is a Symbol and the other either of the allowed types
2282 if isinstance(arg0,Symbol):
2283 sh0=arg0.getShape()
2284 sh1=getShape(arg1)
2285 if not sh0[:axis_offset]==sh1[:axis_offset]:
2286 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))
2287 if isinstance(arg1,float):
2288 return arg0*arg1
2289 elif isinstance(arg1,numpy.ndarray) or isinstance(arg1, Symbol):
2290 return arg0.transposedTensorProduct(arg1, axis_offset)
2291 elif isinstance(arg1, escript.Data):
2292 raise TypeError("tensor product of Symbol and Data not supported yet")
2293 elif isinstance(arg0,numpy.ndarray):
2294 if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
2295 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))
2296 arg0_c=arg0.copy()
2297 arg1_c=arg1.copy()
2298 sh0,sh1=arg0.shape,arg1.shape
2299 d0,d1,d01=1,1,1
2300 for i in sh0[axis_offset:]: d0*=i
2301 for i in sh1[axis_offset:]: d1*=i
2302 for i in sh0[:axis_offset]: d01*=i
2303 arg0_c.resize((d01,d0))
2304 arg1_c.resize((d01,d1))
2305 out=numpy.zeros((d0,d1),numpy.float64)
2306 for i0 in range(d0):
2307 for i1 in range(d1):
2308 out[i0,i1]=numpy.sum(arg0_c[:,i0]*arg1_c[:,i1])
2309 out.resize(sh0[axis_offset:]+sh1[axis_offset:])
2310 return out
2311 elif isinstance(arg0,escript.Data):
2312 if isinstance(arg1, Symbol):
2313 raise TypeError("tensor product of Data and Symbol not supported yet")
2314 # this call has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2315 return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset)
2316
2317 # this should be escript._generalTransposedTensorProduct
2318 def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset):
2319 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2320 return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
2321
2322 def matrix_transposed_mult(arg0,arg1):
2323 """
2324 matrix-transposed(matrix) product of the two arguments.
2325
2326 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[s1,r0]}
2327
2328 The function call ``matrix_transposed_mult(arg0,arg1)`` is equivalent to
2329 ``matrix_mult(arg0,transpose(arg1))``.
2330
2331 The last dimensions of ``arg0`` and ``arg1`` must match.
2332
2333 :param arg0: first argument of rank 2
2334 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2335 :param arg1: second argument of rank 1 or 2
2336 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2337 :return: the product of ``arg0`` and the transposed of ``arg1`` at each data
2338 point
2339 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2340 :raise ValueError: if the shapes of the arguments are not appropriate
2341 """
2342 sh0=getShape(arg0)
2343 sh1=getShape(arg1)
2344 if not len(sh0)==2 :
2345 raise ValueError("first argument must have rank 2")
2346 if not len(sh1)==2 and not len(sh1)==1:
2347 raise ValueError("second argument must have rank 1 or 2")
2348 return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
2349
2350 def tensor_transposed_mult(arg0,arg1):
2351 """
2352 The tensor product of the first and the transpose of the second argument.
2353
2354 For ``arg0`` of rank 2 this is
2355
2356 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[s1,r0]}
2357
2358 and for ``arg0`` of rank 4 this is
2359
2360 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]}
2361
2362 or
2363
2364 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]}
2365
2366 In the first case the the second dimension of ``arg0`` and ``arg1`` must
2367 match and in the second case the two last dimensions of ``arg0`` must match
2368 the two last dimensions of ``arg1``.
2369
2370 The function call ``tensor_transpose_mult(arg0,arg1)`` is equivalent to
2371 ``tensor_mult(arg0,transpose(arg1))``.
2372
2373 :param arg0: first argument of rank 2 or 4
2374 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`
2375 :param arg1: second argument of shape greater of 1 or 2 depending on rank
2376 of ``arg0``
2377 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`
2378 :return: the tensor product of the transposed of ``arg0`` and ``arg1`` at
2379 each data point
2380 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2381 """
2382 sh0=getShape(arg0)
2383 sh1=getShape(arg1)
2384 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2385 return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
2386 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2387 return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
2388 else:
2389 raise ValueError("first argument must have rank 2 or 4")
2390
2391 def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
2392 """
2393 Generalized tensor product of ``arg0`` and transpose of ``arg1``.
2394
2395 C{out[s,t]=Sigma_r arg0[s,r]*arg1[t,r]}
2396
2397 where
2398 - s runs through ``arg0.Shape[:arg0.ndim-axis_offset]``
2399 - r runs through ``arg0.Shape[arg1.ndim-axis_offset:]``
2400 - t runs through ``arg1.Shape[arg1.ndim-axis_offset:]``
2401
2402 The function call ``generalTensorTransposedProduct(arg0,arg1,axis_offset)``
2403 is equivalent to
2404 ``generalTensorProduct(arg0,transpose(arg1,arg1.ndim-axis_offset),axis_offset)``.
2405
2406 :param arg0: first argument
2407 :type arg0: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2408 :param arg1: second argument
2409 :type arg1: ``numpy.ndarray``, `escript.Data`, `Symbol`, ``float``, ``int``
2410 :return: the general tensor product of ``arg0`` and ``transpose(arg1)`` at
2411 each data point
2412 :rtype: ``numpy.ndarray``, `escript.Data`, `Symbol` depending on the input
2413 """
2414 if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2415 arg0,arg1=matchType(arg0,arg1)
2416 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data,
2417 # or one is a Symbol and the other either of the allowed types
2418 if isinstance(arg0,Symbol):
2419 sh0=arg0.getShape()
2420 sh1=getShape(arg1)
2421 r1=getRank(arg1)
2422 if not sh0[arg0.getRank()-axis_offset:]==sh1[r1-axis_offset:]:
2423 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))
2424 if isinstance(arg1,float):
2425 return arg0*arg1
2426 elif isinstance(arg1,numpy.ndarray) or isinstance(arg1, Symbol):
2427 return arg0.tensorTransposedProduct(arg1, axis_offset)
2428 elif isinstance(arg1, escript.Data):
2429 raise TypeError("tensor product of Symbol and Data not supported yet")
2430 elif isinstance(arg0,numpy.ndarray):
2431 if not arg0.shape[arg0.ndim-axis_offset:]==arg1.shape[arg1.ndim-axis_offset:]:
2432 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))
2433 arg0_c=arg0.copy()
2434 arg1_c=arg1.copy()
2435 sh0,sh1=arg0.shape,arg1.shape
2436 d0,d1,d01=1,1,1
2437 for i in sh0[:arg0.ndim-axis_offset]: d0*=i
2438 for i in sh1[:arg1.ndim-axis_offset]: d1*=i
2439 for i in sh1[arg1.ndim-axis_offset:]: d01*=i
2440 arg0_c.resize((d0,d01))
2441 arg1_c.resize((d1,d01))
2442 out=numpy.zeros((d0,d1),numpy.float64)
2443 for i0 in range(d0):
2444 for i1 in range(d1):
2445 out[i0,i1]=numpy.sum(arg0_c[i0,:]*arg1_c[i1,:])
2446 out.resize(sh0[:arg0.ndim-axis_offset]+sh1[:arg1.ndim-axis_offset])
2447 return out
2448 elif isinstance(arg0,escript.Data):
2449 if isinstance(arg1, Symbol):
2450 raise TypeError("tensor product of Data and Symbol not supported yet")
2451 # this call has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2452 return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset)
2453
2454 # this should be escript._generalTensorTransposedProduct
2455 def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset):
2456 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2457 return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
2458
2459 #=========================================================
2460 # functions dealing with spatial dependency
2461 #=========================================================
2462 def grad(arg,where=None):
2463 """
2464 Returns the spatial gradient of ``arg`` at ``where``.
2465
2466 If ``g`` is the returned object, then
2467
2468 - if ``arg`` is rank 0 ``g[s]`` is the derivative of ``arg`` with respect to
2469 the ``s``-th spatial dimension
2470 - if ``arg`` is rank 1 ``g[i,s]`` is the derivative of ``arg[i]`` with
2471 respect to the ``s``-th spatial dimension
2472 - if ``arg`` is rank 2 ``g[i,j,s]`` is the derivative of ``arg[i,j]`` with
2473 respect to the ``s``-th spatial dimension
2474 - if ``arg`` is rank 3 ``g[i,j,k,s]`` is the derivative of ``arg[i,j,k]``
2475 with respect to the ``s``-th spatial dimension.
2476
2477 :param arg: function of which the gradient is to be calculated. Its rank
2478 has to be less than 3.
2479 :type arg: `escript.Data` or `Symbol`
2480 :param where: FunctionSpace in which the gradient is calculated.
2481 If not present or ``None`` an appropriate default is used.
2482 :type where: ``None`` or `escript.FunctionSpace`
2483 :return: gradient of ``arg``
2484 :rtype: `escript.Data` or `Symbol`
2485 """
2486 if isinstance(arg,Symbol):
2487 if where is None:
2488 return arg.grad()
2489 else:
2490 return arg.grad(where)
2491 elif isinstance(arg,escript.Data):
2492 if where==None:
2493 return arg._grad()
2494 else:
2495 return arg._grad(where)
2496 else:
2497 raise TypeError("grad: Unknown argument type.")
2498
2499 def grad_n(arg, n, where=None):
2500 return grad(arg, where)[n]
2501
2502 def integrate(arg,where=None):
2503 """
2504 Returns the integral of the function ``arg`` over its domain. If ``where`` is
2505 present ``arg`` is interpolated to ``where`` before integration.
2506
2507 :param arg: the function which is integrated
2508 :type arg: `escript.Data` or `Symbol`
2509 :param where: FunctionSpace in which the integral is calculated.
2510 If not present or ``None`` an appropriate default is used.
2511 :type where: ``None`` or `escript.FunctionSpace`
2512 :return: integral of ``arg``
2513 :rtype: ``float``, ``numpy.ndarray`` or `Symbol`
2514 """
2515 if isinstance(arg,escript.Data):
2516 if not where==None: arg=escript.Data(arg,where)
2517 if arg.getRank()==0:
2518 return arg._integrateToTuple()[0]
2519 else:
2520 return numpy.array(arg._integrateToTuple())
2521 elif isinstance(arg,Symbol):
2522 return symfn.integrate(arg, where)
2523 else:
2524 arg2=escript.Data(arg,where)
2525 if arg2.getRank()==0:
2526 return arg2._integrateToTuple()[0]
2527 else:
2528 return numpy.array(arg2._integrateToTuple())
2529
2530 def interpolate(arg,where):
2531 """
2532 Interpolates the function into the `FunctionSpace` ``where``. If the
2533 argument ``arg`` has the requested function space ``where`` no interpolation
2534 is performed and ``arg`` is returned.
2535
2536 :param arg: interpolant
2537 :type arg: `escript.Data` or `Symbol`
2538 :param where: `FunctionSpace` to be interpolated to
2539 :type where: `escript.FunctionSpace`
2540 :return: interpolated argument
2541 :rtype: ``escript.Data`` or `Symbol`
2542 """
2543 if isinstance(arg,escript.Data):
2544 if arg.isEmpty():
2545 return arg
2546 elif where == arg.getFunctionSpace():
2547 return arg
2548 else:
2549 return escript.Data(arg,where)
2550 elif isinstance(arg,Symbol):
2551 return symfn.interpolate(arg, where)
2552 else:
2553 return escript.Data(arg,where)
2554
2555 def div(arg,where=None):
2556 """
2557 Returns the divergence of ``arg`` at ``where``.
2558
2559 :param arg: function of which the divergence is to be calculated. Its
2560 shape has to be (d,) where d is the spatial dimension.
2561 :type arg: `escript.Data` or `Symbol`
2562 :param where: `FunctionSpace` in which the divergence will be calculated.
2563 If not present or ``None`` an appropriate default is used.
2564 :type where: ``None`` or `escript.FunctionSpace`
2565 :return: divergence of ``arg``
2566 :rtype: `escript.Data` or `Symbol`
2567 """
2568 if isinstance(arg,escript.Data):
2569 dim=arg.getDomain().getDim()
2570 if not arg.getShape()==(dim,):
2571 raise ValueError("div: expected shape is (%s,)"%dim)
2572 elif not isinstance(arg, Symbol):
2573 raise TypeError("div: argument type not supported")
2574 return trace(grad(arg,where))
2575
2576 def jump(arg,domain=None):
2577 """
2578 Returns the jump of ``arg`` across the continuity of the domain.
2579
2580 :param arg: argument
2581 :type arg: `escript.Data` or `Symbol`
2582 :param domain: the domain where the discontinuity is located. If domain is
2583 not present or equal to ``None`` the domain of ``arg`` is used.
2584 :type domain: ``None`` or `escript.Domain`
2585 :return: jump of ``arg``
2586 :rtype: `escript.Data` or `Symbol`
2587 """
2588 if domain==None: domain=arg.getDomain()
2589 return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
2590
2591 def L2(arg):
2592 """
2593 Returns the L2 norm of ``arg`` at ``where``.
2594
2595 :param arg: function of which the L2 norm is to be calculated
2596 :type arg: `escript.Data` or `Symbol`
2597 :return: L2 norm of ``arg``
2598 :rtype: `float` or `Symbol`
2599 :note: L2(arg) is equivalent to ``sqrt(integrate(inner(arg,arg)))``
2600 """
2601 if isinstance(arg,Symbol):
2602 return symfn.L2(arg)
2603 return sqrt(integrate(inner(arg,arg)))
2604
2605 def getClosestValue(arg,origin=0):
2606 """
2607 Returns the value in ``arg`` which is closest to origin.
2608
2609 :param arg: function
2610 :type arg: `escript.Data`
2611 :param origin: reference value
2612 :type origin: ``float`` or `escript.Data`
2613 :return: value in ``arg`` closest to origin
2614 :rtype: ``numpy.ndarray``
2615 """
2616 return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
2617
2618 def normalize(arg,zerolength=0):
2619 """
2620 Returns the normalized version of ``arg`` (=``arg/length(arg)``).
2621
2622 :param arg: function
2623 :type arg: `escript.Data` or `Symbol`
2624 :param zerolength: relative tolerance for arg == 0
2625 :type zerolength: ``float``
2626 :return: normalized ``arg`` where ``arg`` is non-zero, and zero elsewhere
2627 :rtype: `escript.Data` or `Symbol`
2628 """
2629 l=length(arg)
2630 m=whereZero(l,zerolength*Lsup(l))
2631 mm=1-m
2632 return arg*(mm/(l*mm+m))
2633
2634 def deviatoric(arg):
2635 """
2636 Returns the deviatoric version of ``arg``.
2637 """
2638 return arg-(trace(arg)/trace(kronecker(arg.getDomain())))*kronecker(arg.getDomain())
2639
2640 def vol(arg):
2641 """
2642 Returns the volume or area of the oject ``arg``
2643
2644 :param arg: a geometrical object
2645 :type arg: `escript.FunctionSpace` or `escript.Domain`
2646 :rtype: ``float``
2647 """
2648 if isinstance(arg,escript.Domain): arg=escript.Function(arg)
2649 return integrate(escript.Scalar(1.,arg))
2650
2651 def meanValue(arg):
2652 """
2653 return the mean value of the argument over its domain
2654
2655 :param arg: function
2656 :type arg: `escript.Data`
2657 :return: mean value
2658 :rtype: ``float`` or {numpy.ndarray}
2659 """
2660 fs=arg.getFunctionSpace()
2661 d=fs.getDomain()
2662 if fs == escript.Solution(d) or fs == escript.ContinuousFunction(d):
2663 fs=escript.Function(d)
2664 if fs == escript.ReducedSolution(d) or fs == escript.ReducedContinuousFunction(d):
2665 fs=escript.ReducedFunction(d)
2666 a=vol(fs)
2667 if a == 0:
2668 raise ValueError("FunctionSpace %s with zero volume."%str(fs))
2669 return integrate(arg,fs)/a
2670
2671 def diameter(domain):
2672 """
2673 Returns the diameter of a domain.
2674
2675 :param domain: a domain
2676 :type domain: `escript.Domain`
2677 :rtype: ``float``
2678 """
2679 out=0
2680 for v in boundingBox(domain): out+=(v[1]-v[0])**2
2681 return sqrt(out)
2682
2683 def boundingBox(domain):
2684 """
2685 Returns the bounding box of a domain
2686
2687 :param domain: a domain
2688 :type domain: `escript.Domain`
2689 :return: bounding box of the domain
2690 :rtype: ``list`` of pairs of ``float``
2691 """
2692 x=domain.getX()
2693 out=[]
2694 for i in range(domain.getDim()):
2695 x_i=x[i]
2696 out.append((inf(x_i),sup(x_i)))
2697 return out
2698
2699 def longestEdge(domain):
2700 """
2701 Returns the length of the longest edge of the domain
2702
2703 :param domain: a domain
2704 :type domain: `escript.Domain`
2705 :return: longest edge of the domain parallel to the Cartisean axis
2706 :rtype: ``float``
2707 """
2708 return max([v[1]-v[0] for v in boundingBox(domain) ])
2709
2710 def mkDir(*pathname):
2711 """
2712 creates a directory of name ``pathname`` if the directory does not exist.
2713
2714 :param pathname: valid path name
2715 :type pathname: ``str`` or ``sequence of strings``
2716 :note: The method is MPI safe.
2717 """
2718 errno=0
2719 p_fail=None
2720 if getMPIRankWorld()==0:
2721 for p in pathname:
2722 if os.path.exists(p):
2723 if not os.path.isdir(p):
2724 errno=2
2725 p_fail=p
2726 else:
2727 try:
2728 os.makedirs(p)
2729 except Exception as e:
2730 errno=1
2731 p_fail=p
2732
2733 errno=getMPIWorldMax(errno)
2734 if errno>0:
2735 if errno==2:
2736 if p_fail == None:
2737 raise IOError("Unable to create directory.")
2738 else:
2739 raise IOError("Unable to create directory %s. It already exists and is not a directory."%p_fail)
2740 elif e==None:
2741 if p_fail == None:
2742 raise IOError("Unable to create directory.")
2743 else:
2744 raise IOError("Unable to create directory %s."%p_fail)
2745 else:
2746 if hasattr(e,"message"):
2747 raise IOError(e.message)
2748 else:
2749 if p_fail == None:
2750 raise IOError("Unable to create directory.")
2751 else:
2752 raise IOError("Unable to create directory %s."%p_fail)
2753
2754 class FileWriter(object):
2755 """
2756 Interface to write data to a file. In essence this class wrappes the standard ``file`` object to write data that are global in MPI
2757 to a file. In fact, data are writen on the processor with MPI rank 0 only. It is recommended to use ``FileWriter`` rather than ``open`` in order to write
2758 code that is running with as well as with MPI. It is safe to use ``open`` onder MPI to read data which are global under MPI.
2759 :var name: name of file
2760 :var mode: access mode (='w' or ='a')
2761 :var closed: True to indicate closed file
2762 :var newlines: line seperator
2763 """
2764 def __init__(self,fn,append=False,createLocalFiles=False):
2765 """
2766 Opens a file of name ``fn`` for writing. If running under MPI only the first processor with rank==0
2767 will open the file and write to it. If ``createLocalFiles`` each individual processor will create a file
2768 where for any processor with rank>0 the file name is extended by its rank. This option is normally only used for
2769 debug purposes.
2770
2771 :param fn: filename.
2772 :type fn: ``str``
2773 :param append: switches on the creation of local files.
2774 :type append: ``bool``
2775 :param createLocalFiles: switches on the creation of local files.
2776 :type createLocalFiles: ``bool``
2777 """
2778 errno=0
2779 self.name=fn
2780 if append:
2781 self.mode='a'
2782 else:
2783 self.mode='w'
2784 self.__file=None
2785 self.closed=False
2786 self.newlines=os.linesep
2787 e=None
2788 # if not the master:
2789 if getMPIRankWorld()>0:
2790 if createLocalFiles:
2791 fn2=fn+".%s"%getMPIRankWorld()
2792 try:
2793 self.__file=open(fn2,self.mode)
2794 except Exception as e:
2795 errno=1
2796 else:
2797 try:
2798 self.__file=open(fn,self.mode)
2799 except Exception as e:
2800 errno=1
2801 self.__handelerror(errno,e,"opening")
2802
2803 def __handelerror(self,errno,e,operation):
2804 errno=getMPIWorldMax(errno)
2805 if errno>0:
2806 if e==None:
2807 raise IOError("Unable to access file %s in mode %s for %s."%(self.name,self.mode,operation))
2808 else:
2809 raise IOError(str(e))
2810
2811 def close(self):
2812 """
2813 Closes the file
2814 """
2815 errno=0
2816 e=None
2817 try:
2818 if not self.__file == None:
2819 self.__file.close()
2820 except Exception as e:
2821 errno=1
2822 self.__handelerror(errno,e,"closing")
2823 self.closed=True
2824
2825 def flush(self):
2826 """
2827 Flush the internal I/O buffer.
2828 """
2829 errno=0
2830 e=None
2831 try:
2832 if not self.__file == None:
2833 self.__file.flush()
2834 except Exception as e:
2835 errno=1
2836 self.__handelerror(errno,e,"flushing")
2837
2838 def write(self,txt):
2839 """
2840 Write string ``txt`` to file.
2841
2842 :param txt: string ``txt`` to be written to file
2843 :type txt: ``str``
2844 """
2845 errno=0
2846 e=None
2847 try:
2848 if not self.__file == None:
2849 self.__file.write(txt)
2850 except Exception as e:
2851 errno=1
2852 self.__handelerror(errno,e,"writing")
2853
2854 def writelines(self, txts):
2855 """
2856 Write the list ``txt`` of strings to the file.
2857
2858 :param txts: sequense of strings to be written to file
2859 :type txts: any iterable object producing strings
2860 :note: Note that newlines are not added. This method is equivalent to call write() for each string.
2861 """
2862 errno=0
2863 e=None
2864 try:
2865 if not self.__file == None:
2866 self.__file.writelines(txts)
2867 except Exception as e:
2868 errno=1
2869 self.__handelerror(errno,e,"writing strings")
2870
2871 def reorderComponents(arg,index):
2872 """
2873 Resorts the components of ``arg`` according to index.
2874
2875 """
2876 raise NotImplementedError
2877
2878 def showEscriptParams():
2879 """
2880 Displays the parameters escript recognises with an explanation and their
2881 current value.
2882 """
2883 p=listEscriptParams()
2884 for name,value,desc in p:
2885 print('%s (=%s): %s'%(name, value, desc))
2886
2887 #Lazy related things
2888 #These are just wrappers
2889 def resolve(arg):
2890 """
2891 Returns the value of arg resolved.
2892 """
2893 if not isinstance(arg,Data):
2894 raise TypeError("Can only resolve Data.")
2895 if arg.isLazy():
2896 arg.resolve()
2897 return arg
2898
2899 def delay(arg):
2900 """
2901 Returns a lazy version of arg
2902 """
2903 if not isinstance(arg,Data):
2904 raise TypeError("Can only delay Data.")
2905 return arg.delay()
2906
2907 def positive(arg):
2908 """
2909 returns the positive part of arg
2910 """
2911 return (abs(arg)+arg)/2.
2912
2913 def negative(arg):
2914 """
2915 returns the negative part of arg
2916 """
2917 return (arg-abs(arg))/2.
2918
2919 def safeDiv(arg0, arg1, rtol=None):
2920 """
2921 returns arg0/arg1 but return 0 where arg1 is (almost) zero
2922 """
2923 if rtol==None:
2924 m1=whereZero(arg1,tol=0)
2925 else:
2926 m1=whereZero(arg1,tol=None, rtol=rtol)
2927 return arg0/(arg1+m1)*(1-m1)
2928
2929 def condEval(f, tval, fval):
2930 """
2931 Wrapper to allow non-data objects to be used.
2932 """
2933 if not isinstance(tval,Data) and not isinstance(fval,Data):
2934 raise TypeError("At least one of the alternatives must be a Data object.")
2935 if isinstance(tval,Data) and isinstance(fval, Data):
2936 return _condEval(f,tval,fval)
2937 if not isinstance(fval, Data):
2938 return _condEval(f, tval, Data(fval, tval.getShape(), tval.getFunctionSpace()))
2939 return _condEval(f, Data(fval, fval.getShape(), fval.getFunctionSpace()), fval )
2940

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26