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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26