/[escript]/trunk/escript/py_src/util.py
ViewVC logotype

Contents of /trunk/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26