Source code for geotecha.mathematics.multi_transform
# geotecha - A software suite for geotechncial engineering
# Copyright (C) 2018 Rohan T. Walker (rtrwalker@gmail.com)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see http://www.gnu.org/licenses/gpl.html.
"""Some routines related multi dimensional integral transforms"""
from __future__ import division, print_function
import numpy as np
from matplotlib import pyplot as plt
import functools
import unittest
from numpy.testing import assert_allclose
from geotecha.mathematics.laplace import Talbot
from geotecha.mathematics.hankel import HankelTransform
from geotecha.mathematics.fourier import FourierTransform
import collections
[docs]def ntransform(func, transforms, transvars, args=None, opts=None):
"""
Multi-dimensional integral transforms over multiple variables.
General idea and organisation of code is from scipy.integrate.nquad
Parameters
----------
func : callable
The function to be integrated. Has arguments of ``x0, ... xn``,
``t0, tm``, where integration is carried out over ``x0, ... xn``, which
must be floats. Function signature should be
``func(x0, x1, ..., xn, t0, t1, ..., tm)``. Integral transforms are
carried out in order. That is, integration over ``x0`` is the
innermost integral/transform, and ``xn`` is the outermost.
transform : iterable object of string
Each elmement of transforms is a string corresponding to one of the
available transforms. Current options are:
- Hankel
- Hankel_inverse
- Fourier
- Fourier_inverse
- Laplace_inverse
transvars : iterable
Transformation variable for each transformation.
args : iterable object, optional
Additional arguments ``t0, ..., tn``, required by `func`.
opts : iterable object or dict, optional
Options to be passed to each transform. May be empty, a dict, or
a sequence of dicts, or functions that return a dict. If empty, the
default options of each transform are used. If a dict, the same
options are used for all levels of integration. If a sequence, then each
element of the sequence corresponds to a particular integration/
transform. e.g.
opts[0] corresponds to integration/transform over x0, and so on.
See the individual transforms for options.
Returns
-------
result : float
The result of the integration.
abserr : float
The maximum of the estimates of the absolute error in the various
integration results.
See Also
--------
geotecha.mathematics.laplace.Talbot : numerical inverse Laplace
geotecha.mathematics.hankel.HankelTransform : Hankel transform
geotecha.mathematics.fourier.FourierTransform : Fourier transform
Notes
-----
ntransform is quite temperamental. Be careful with the following:
- For some reason performing repeated multi dimensional fourier
transforms with the default integration limit of np.inf does not
work (I get a very small number... the wrong number). This problem
also happens when doing the same transforms using scipy's nquad for
multidemensional integration. There is some issue with the recursion
and use of the underlying QUADPACK integration routines when performing
integrations with a cos or sin weight function with infinte integration
limits. To get around this you could try truncating the integral by
setting the `b` keyword agument to the FourierTransform object. You
will probably have to trial a few values; larger values give more
accuracy but too large and the solution is gibberish.
- Because the inverse laplace transform uses imaginary and any QUADPACK
integations can only use real numbers. A laplace transfrom usually must
be the inner most intergal, i.e. first in the `transform` list.
If performing multi-dimensional inverse lapace transforms then
they must be at the front of the `transform` list; all inverse
laplace transforms other than the first in the `transform` list must
have {'vectorised': False}, in the corresponding `opts` dict.
If `f` is not vectorised then the first inverse Laplace transform must
also have {'vectorised': False}
- Because the Hankel transform uses numpy broadcasting, it must always
be the inner most integral, i.e. first in the `transform` list. If also
doing an inverse laplace transform, then 'Laplace_inverse` must be
second in the `transform` list with {'vectorised': False} in the
corresponding `opts` dict.
- Rules of thumb are 1. 'Hankel_transform', if it occurs, must be first
in the `transform` list, 2. 'Laplace_inverse' must be before any
'Fourier' or 'Fourier_inverse' instance in the `transform` list, 3. For
any but the first 'Laplace_inverse' instance in the `transform` list,
the corresponding `opts` dict must contain 'vectorized':False, 4.
'Fourier' or 'Fourier_inverse' may fail using the default b=np.inf, try
truncating the integral with 'b':35 in the corresponding `opts` dict.
- Make sure that the order of the `transform` list matches the order
or the args in function `f`. First in list transforms first arg of
`f`, 2nd in list transforms secnod arg of `f`.
"""
depth = len(transforms)
try:
assert len(transvars)==depth
except TypeError:
raise TypeError('transvar must be iterable')
except AssertionError:
raise AssertionError('transvar and transforms must be the same length')
if args is None:
args = ()
if opts is None:
opts = [dict([])] * depth
if isinstance(opts, dict):
opts = [opts] * depth
else:
opts = [opt if isinstance(opt, collections.Callable) else
_OptFunc(opt) for opt in opts]
return _NTransform(func, transforms,
transvars, opts).integral_transform(*args)
class _OptFunc(object):
"""When called the object will return the variable used to initialize it
Helper function for the ntransform function. Useful when you are using a
function that only accepts callable objects. Say your function expects
a callable object that returns a two element list. If you always want
to provide a known two element list then you could hard code a
function to return your list, or you could pass _OptFunc(my_list) and
avoid all the code snippets.
Parameters
----------
opt : anything
The variable that will be returned when the object is called with
any positional arguments.
"""
def __init__(self, opt):
self.opt = opt
def __call__(self, *args):
"""Return stored dict."""
return self.opt
class _NTransform(object):
"""Recursively perform integral transforms
Engine room of the ntransform function.
"""
def __init__(self, func, transforms, transvars, opts):
self.abserr = 0
self.func = func
self.transforms = transforms
self.transvars = transvars
self.opts = opts
self.maxdepth = len(transforms)
self.tdict= {'Hankel': self.hankel,
'Hankel_inverse': self.hankel,
'Fourier': self.fourier,
'Fourier_inverse': self.fourier_inverse,
'Laplace_inverse': self.laplace_inverse,}
def hankel(self, f, tvar, args, **opt):
ht = HankelTransform(f, args, **opt)
val, err = ht(tvar)
return val, err
def fourier(self, f, tvar, args, **opt):
ft = FourierTransform(f, args, **opt)
val, err = ft(tvar)
return val, err
def fourier_inverse(self,f, tvar, args, **opt):
ft = FourierTransform(f, args, inv=True, **opt)
val, err = ft(tvar)
return val, err
def laplace_inverse(self, f, tvar, args, **opt):
ilt = Talbot(f, **opt)
val = ilt(tvar, args)
return val, 0.0
def integral_transform(self, *args, **kwargs):
"""Perform the transforms"""
depth = kwargs.pop('depth', 0)
if kwargs:
raise ValueError('unexpected kwargs')
# Get the integration range and options for this depth.
ind = -(depth + 1)
transform = self.tdict[self.transforms[ind]]
tvar = self.transvars[ind]
fn_opt = self.opts[ind]
opt = dict(fn_opt(*args))
if depth + 1 == self.maxdepth:
f = self.func
else:
f = functools.partial(self.integral_transform, depth=depth+1)
value, abserr = transform(f, tvar, args=args, **opt)
self.abserr = max(self.abserr, np.abs(abserr))
if depth > 0:
return value
else:
# Final result of n-D integration with error
return value, self.abserr
#def f_f1(x, y, a, b):
# f1 = np.exp(-a * abs(x))
# f2 = np.exp(-b * abs(y))
# return f1*f2
#def f_f1_(x, y, a, b):
# f1 = 2 * a / (a**2 + x**2)
# f2 = 2 * b / (b**2 + y**2)
# return f1*f2
#
#def test_f_f1():
# a, b, c = 1.8, 2.2, 0.5
# x, y, z = 1.1, 1.2, 1.3
## a, b, c = 0.15, 0.18,0.22
## x, y, z = 0.1, 0.2, 0.3
# tvar = (x, y)
# args= (a, b)
# transforms=['Fourier', 'Fourier']
# opts=[{'func_is_real': True, 'real_part_even': True, 'b': 10}]*2
# print(f_f1_(*(tvar+args)))
# print(f_f1_(x, y, a, b))
# assert_allclose(ntransform(f_f1, transforms, tvar, args, opts)[0],
# f_f1_(*(tvar+args)), atol=0)
#
#
#
#
#def aaa():
# a, b= 0.15, 0.18
# x, y = 0.1, 0.2
# a, b= 1, 1
# x, y = 20, 1
# from scipy.integrate import nquad
# from scipy.integrate import quad
#
# tvar = (x, y)
# args= (a, b)
# transforms=['Fourier', 'Fourier']
# opts=[{'func_is_real': True, 'real_part_even':True, 'b':10}]*2
# val, err=ntransform(f_f1, transforms, tvar, args, opts)
# expected = f_f1_(*(tvar+args))
# print('expected', expected)
# print('ntransform',val, err)
#
# epsrel=1.49e-08
# epsabs=1.49e-17
# val, err = nquad(f_f1, [(0, np.inf), (0, 54)], args=(a, b),
# opts=[{'weight': 'cos', 'wvar': x, 'epsabs':epsabs, 'epsrel': epsrel}
# , {'weight': 'cos', 'wvar': y, 'epsabs':epsabs, 'epsrel': epsrel}])
# print('nquad', 4*val, err)
# val,err = quad(lambda x,a: np.exp(-a * abs(x)) , 0, np.inf, args=(a,),
# weight='cos', wvar=x)
# val1, err1 = quad(lambda x,a: np.exp(-a * abs(x)) , 0, np.inf, args=(b,),
# weight='cos', wvar=y)
# print('quad*quad', 4*val*val1, err+err1)
if __name__ == '__main__':
# aaa()
# test_f_f1()
import nose
nose.runmodule(argv=['nose', '--verbosity=3', '--with-doctest'])
# nose.runmodule(argv=['nose', '--verbosity=3'])