Source code for geotecha.mathematics.fourier

# 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.
"""Fourier transforms."""


from __future__ import division, print_function
import matplotlib.pyplot
import numpy as np
from scipy import integrate
from scipy.special import jn_zeros
from scipy.special import jn
from matplotlib import pyplot as plt
import functools

import unittest
from numpy.testing import assert_allclose


from geotecha.mathematics.quadrature import gl_quad
from geotecha.mathematics.quadrature import gk_quad
from geotecha.mathematics.quadrature import gauss_kronrod_abscissae_and_weights
from geotecha.mathematics.quadrature import gauss_legendre_abscissae_and_weights
from geotecha.mathematics.quadrature import shanks_table
from geotecha.mathematics.quadrature import shanks


[docs]def real_func(x, *myargs): """Real part of a function Basically return np.real(func(x, *myargs[1:])) where func is the first argument after x. Parameters ---------- x : float Value to evaluate function at. func : function/callable Function from which to return the real part. Always the first argument after `x`. myargs : optional Any remaining arguments will be passed to func(x, *myargs[1:]). Returns ------- out : ndarray Real part of func(x, *myargs) See Also -------- imag_func : imaginary part of function Examples -------- >>> def f(x, a): ... return a*x+a*1.j >>> real_func(2, f, 4) 8.0 >>> real_func(3.j,f, 2) 0.0 >>> real_func(np.array([3.j, 1+2.j]),f, 2) array([0., 2.]) """ func = myargs[0] myargs = (x,) + myargs[1:] return +np.real(func(*myargs))
[docs]def imag_func(x, *myargs): """Imaginary part of a function Basically return np.imag(func(x, *myargs[1:])) where func is the first argument after x. Parameters ---------- x : float Value to evaluate function at. func : function/callable Function from which to return the imaginary part. Always the first argument after `x` myargs : optional Any remaining arguments will be passed to func(x, *myargs[1:]). Returns ------- out : ndarray Imaginary part of func(x, *myargs). See Also -------- real_func : real part of function Examples -------- >>> def f(x, a): ... return a*x+a*1.j >>> imag_func(2, f, 4) 4.0 >>> imag_func(3.j,f, 2) 8.0 >>> imag_func(np.array([3.j, 2.j]),f, 2) array([8., 6.]) """ func = myargs[0] myargs = (x,) + myargs[1:] return +np.imag(func(*myargs))
[docs]def func_mirror_for_even_weight(x, *myargs): """Mirror a function abount the y-axis Given a composite function f(x) * w(x) where w(x) is an even weighting function, return g(x) such that g(x)*w(x) gives same value as f(-x)*w(-x). This can be useful in transforming a fourier cosine integral with negative integation limits to one with positive limits. Parameters ---------- x : float Value to evaluate function at. func : function/callable Function to mirror. Always the first argument after `x`. myargs : optional Any remaining arguments will be passed to `func`. Returns ------- out : ndarray Value of func(-x, *myargs). See Also -------- func_mirror_for_odd_weight : mirror for an odd weight function Examples -------- >>> def f(x, a): ... return a*x+1 >>> func_mirror_for_even_weight(5, f, 2) -9 >>> def ff(x, a): ... return a*x + 1.j >>> func_mirror_for_even_weight(3, real_func, ff, 4) -12.0 """ y = -x func = myargs[0] myargs = (y,) + myargs[1:] return +func(*myargs)
[docs]def func_mirror_for_odd_weight(x, *myargs): """Rotate function by 180 degrees (or mirror about x and y-axis in turn). Given a composite function f(x) * w(x) where w(x) is an odd weighting function, return g(x) such that g(x)*w(x) gives same value as f(-x)*w(-x). This can be useful in transforming a fourier sine integral with negative integration limits to one with positive limits. Parameters ---------- x : float Value to evaluate function at. func : function/callable Function to mirror. Always the first argument after `x`. myargs : optional Any remaining arguments will be passed to `func`. Returns ------- out : ndarray Value of -func(-x, *myargs) See Also -------- func_mirror_for_even_weight : mirror for an even wieght function Examples -------- >>> def f(x, a): ... return a*x+1 >>> func_mirror_for_odd_weight(5, f, 2) 9 >>> def ff(x, a): ... return a*x + 1.j >>> func_mirror_for_odd_weight(3, real_func, ff, 4) 12.0 """ y = -x func = myargs[0] myargs = (y,) + myargs[1:] return -func(*myargs)
[docs]def cosine_transform(func, w, args=(), a=0.0, b=np.inf): """Fourier cosine transform Note that any function that can divide by zero may cause problems because QUADPACK includes the end points in integration. Parameters ---------- func : function/callable Function to transform. `func` will be called func(x, *args). `func` must return a real. w : float Transform variable. args : tuple, optional Arguments to pass to `func` a, b : float, optional Integration limits. Defualt a=0.0, b=np.inf. Returns ------- value : float Value of transform at w err : float Error estimate from quadpack Notes ----- The fourier cosine transform is given by: .. math:: F_c=\mathcal{F}_c\\{f(x)\\}(w) = \\int_0^{\\infty}f(x)\\cos(wx)\\,\\mathrm{d}x """ return integrate.quad(func, a, b, args=args, weight='cos', wvar=w)
[docs]def sine_transform(func, w, args=(), a=0.0, b=np.inf): """Fourier sine transform Note that any function that can divide by zero may cause problems because QUADPACK includes the end points in integration. Parameters ---------- func : function/callable Function to transform. `func` will be called func(x, *args). `func` must return a real. w : float Transform varibale. args : tuple, optional Arguments to pass to `func` a, b : float, optional Integration limits. Defualt a=0.0, b=np.inf. Returns ------- value : float Value of transform at w. err : float Error estimate from quadpack. Notes ----- The fourier sine transform is given by: .. math:: F_s=\mathcal{F}_s\\{f(x)\\}(w) = \\int_0^{\\infty}f(x)\\sin(wx)\\,\\mathrm{d}x """ return integrate.quad(func, a, b, args=args, weight='sin', wvar=w)
[docs]class FourierTransform(object): """One dimensional Fourier transform using scipy.quad Note that any function that can divide by zero may cause problems because QUADPACK includes the end points in integration. Parameters ---------- func : function Function to transform. `func` is called by func(x, *args) args : tuple, optional tuple of arguments to pass to func. Default args=(). inv : True/False, optional If True then the inverse Fourier transform will be performed. Default inv=False. func_is_real : True/False, optional If True then func is purely real. It returns a real number. default func_is_real=False. func_is_imag : True/False, optional If True then func is purely imaginary. It returns a real number that should be multiplied by i. Default func_is_imag=False real_part_even : True/False, optional If True then the real part of func is even. Default real_part_even=False. real_part_odd : True/False, optional If True then the real part of func is odd. Default real_part_odd=False. imag_part_even : True/False, optional If True then the imaginary part of func is even. Default imag_part_even=False. imag_part_odd : True/False, optional If True then the imaginary part of func is odd. Default imag_part_odd=False. a, b : float, optional Integration limits. Defualt a=0.0, b=np.inf Attributes ---------- inv_sign : float The sign of some expressions change for the inverse fourier transform, `inv_sign` accounts for that sign change. If inv=True, inv_sign=-1; If inv=False, inv_sign=+1. inv_const : float For inverse fourier transform all expressions are multiplied by 1/(2*pi). """ def __init__(self, func, args=(), inv=False, func_is_real=False, func_is_imag=False, real_part_even=False, real_part_odd=False, imag_part_even=False, imag_part_odd=False, a=0.0, b=np.inf): self.func = func self.args = args self.inv = inv if self.inv: self.inv_const=(2.0*np.pi)**(-1) #fourier inverse is 1/(2*pi) * integral self.inv_sign = -1 # some fourier inverse terms are negative else: self.inv_const = 1 self.inv_sign = 1 self.func_is_real = func_is_real self.func_is_imag = func_is_imag self.real_part_even = real_part_even self.real_part_odd = real_part_odd self.imag_part_even = imag_part_even self.imag_part_odd = imag_part_odd self.a = a self.b = b self.fargs=(self.func,) + self.args
[docs] def real_func(self, x): """Real part of func""" return +np.real(self.func(x, *self.args))
[docs] def imag_func(self, x): """Imaginary part of func""" return +np.imag(self.func(x, *self.args))
[docs] def mfe(self, x): """Mirror func for even weight function""" return func_mirror_for_even_weight(x, *self.fargs)
[docs] def mfo(self, x): """Mirror func for odd weight function""" return func_mirror_for_odd_weight(x, *self.fargs)
[docs] def mfre(self, x): """Mirror real(func) for even weight function""" return +np.real(func_mirror_for_even_weight(x, *self.fargs))
[docs] def mfro(self, x): """Mirror real(func) for odd weight function""" return +np.real(func_mirror_for_odd_weight(x, *self.fargs))
[docs] def mfie(self, x): """Mirror imag(func) for even weight function""" return +np.imag(func_mirror_for_even_weight(x, *self.fargs))
[docs] def mfio(self, x): """Mirror imag(func) for odd weight function""" return +np.imag(func_mirror_for_odd_weight(x, *self.fargs))
def __call__(self, s): """Perform 1d Fourier transform at s Parameters ---------- s : float Transform variable """ if self.func_is_real: if self.real_part_even: igral, err = [2 * self.inv_const * v for v in cosine_transform(self.func, s, self.args, a=self.a, b=self.b)] return igral, err elif self.real_part_odd: igral, err = [2.j * self.inv_const * self.inv_sign * v for v in sine_transform(self.func, s, self.args, a=self.a, b=self.b)] return -igral, err else: #func is real and exhibits no symmetry. igral = 0 err = 0 # real transform result [-inf, 0] ig, er = [self.inv_const * v for v in cosine_transform(self.mfe, s, a=self.a, b=self.b)] igral += ig; err += er # real transform result [-, inf] ig, er = [self.inv_const * v for v in cosine_transform(self.func, s, self.args, a=self.a, b=self.b)] igral += ig; err += er # imag transform result [0, inf] ig, er = [1.j * self.inv_const * self.inv_sign * v for v in sine_transform(self.mfo, s, a=self.a, b=self.b)] igral -= ig; err += er # imag transform result [-inf, 0] ig, er = [1.j * self.inv_const * self.inv_sign * v for v in sine_transform(self.func, s, self.args, a=self.a, b=self.b)] igral -= ig; err += er return igral, err if self.func_is_imag: if self.imag_part_even: igral, err = [2.j * self.inv_const * v for v in cosine_transform(self.func, s, self.args, a=self.a, b=self.b)] return igral, err elif self.imag_part_odd: igral, err = [2 * self.inv_const * self.inv_sign * v for v in sine_transform(self.func, s, self.args, a=self.a, b=self.b)] return igral, err else: #func is imaginary and ehibits non symmetry igral = 0 err = 0 # imag transform result [-inf, 0] ig, er = [1.j * self.inv_const * v for v in cosine_transform(self.mfe, s, a=self.a, b=self.b)] igral += ig; err += er # imag transform result [0, inf] ig,er = [1.j * self.inv_const * v for v in cosine_transform(self.func, s, self.args, a=self.a, b=self.b)] igral += ig; err += er # real transform result [-inf, 0] ig, er = [self.inv_const * self.inv_sign * v for v in sine_transform(self.mfo, s, a=self.a, b=self.b)] igral -= ig; err += er # real transform result [0, inf] ig, er = [self.inv_const * self.inv_sign * v for v in sine_transform(self.func, s, self.args, a=self.a, b=self.b)] igral -= ig; err += er return igral, err #if we have reached here then func is complex #use real and imag parts of func igral = 0 err = 0 if self.real_part_even: ig, er = [2 * self.inv_const * v for v in cosine_transform(self.real_func, s, a=self.a, b=self.b)] igral += ig; err += er elif self.real_part_odd: # 2 * I(real_func * sin(s*x), 0, +inf) ig, er = [2.j * self.inv_const * self.inv_sign * v for v in sine_transform(self.real_func, s, a=self.a, b=self.b)] igral -= ig; err += er else: #real part of function exhibits no symmetry # real transform result [-inf, 0] ig, er = [self.inv_const * v for v in cosine_transform(self.mfre, s, a=self.a, b=self.b)] igral += ig; err += er # real transform result [0, inf] ig, er = [self.inv_const * v for v in cosine_transform(self.real_func, s, a=self.a, b=self.b)] igral += ig; err += er # imag transform result [-inf, 0] ig, er = [1.j * self.inv_const * self.inv_sign * v for v in sine_transform(self.mfro, s, a=self.a, b=self.b)] igral -= ig; err += er # imag transform result [0, inf] ig, er = [1.j * self.inv_const * self.inv_sign * v for v in sine_transform(self.real_func, s, a=self.a, b=self.b)] igral -= ig; err += er if self.imag_part_even: igral, err = [2.j * self.inv_const * v for v in cosine_transform(self.func, s, self.args, a=self.a, b=self.b)] igral += ig; err += er elif self.real_part_odd: igral, err = [2 * self.inv_const * self.inv_sign * v for v in sine_transform(self.func, s, self.args, a=self.a, b=self.b)] igral += ig; err += er else: #imag part of function exhibits no symmetry # imag transform result [-inf, 0] ig, er = [1.j * self.inv_const * v for v in cosine_transform(self.mfie, s, a=self.a, b=self.b)] igral += ig; err += er # imag transform result [0, inf] ig, er = [1.j * self.inv_const * v for v in cosine_transform(self.imag_func, s, a=self.a, b=self.b)] igral += ig; err += er # real transform result [-inf, 0] ig, er = [self.inv_const * self.inv_sign * v for v in sine_transform(self.mfio, s, a=self.a, b=self.b)] igral += ig; err += er # real transform result [0, inf] ig, er = [self.inv_const * self.inv_sign * v for v in sine_transform(self.imag_func, s, a=self.a, b=self.b)] igral += ig; err += er return igral, err
[docs]def vcosine_transform(f, s, args=(), m=20, ng=20, shanks_ind=None): """Cosine transform of f(x) at transform variable s This is a vectorized cosine transform. Parameters ---------- f : function or method Function to apply cosine trasnform to. f is called with f(x, *args). s : 1d array Coordinate(s) to evaluate transform at. args : tuple, optional arguments to pass to f m : int, optional Number of segments to break the integration interval into. Each segment will be between the zeros of the cos function, Default m=20. ng : [2-20, 32, 64, 100], optional Number of gauss points to use in integration., Default ng=20. shanks_ind : int, optional Start position of intervals to start shanks extrapolation. Default shanks_ind=None i.e. no extrapolation. Be careful when using shanks extrapolation; make sure you only begin to use it after the intgrand is well behaved. Returns ------- f : 1d array of float Value of transform at s Notes ----- Careful with singularities. Because there is no way to increase the integration points at a particular spot the infinite behaviur may not be captured well. For example x**-0.5 should transform to sqrt(pi/2*w) but due to the sinularity at x=0 it does not converge well even using ng=100. """ si = np.atleast_1d(s) xk_, wk = gauss_legendre_abscissae_and_weights(ng) # integration intervals zeros = np.zeros(m + 1, dtype=float) zeros[1:] = (2 * np.arange(m) + 1) * np.pi/2 aj = zeros[0:-1] bj = zeros[1:] #dims of array: # 0 or i dim is each transform coordinate # 1 or j dim is each integration interval # 2 or k dim is each integration point # 2 dim will be summed to get integral of each interval # 1 dim will be summed or shanks'ed to give transform at each coord # si = si[:, np.newaxis, np.newaxis] aj = aj[np.newaxis, :, np.newaxis] bj = bj[np.newaxis, :, np.newaxis] xk_ = xk_[np.newaxis, np.newaxis, :] wk = wk[np.newaxis, np.newaxis, :] aij = aj / si bij = bj / si bma = (bij - aij) / 2 # b minus a bpa = (aij + bij) /2 # b plus a xijk = bma * xk_ + bpa # xj_ are in [-1, 1] so need to transform to [a, b] fijk = f(xijk, *args) fijk *= np.cos(si * xijk) # fijk *= xijk igral = bma[:,:,0] * np.sum(fijk * wk, axis=2) if shanks_ind is None: return igral.sum(axis=1) else: #extrapolate igral.cumsum(axis=1 , out=igral) return shanks(igral, shanks_ind)
[docs]def v2dcosine_transform(f, s1, s2, args=(), m=20, ng=20, shanks_ind=None): """Cosine transform of f(x, y) at transform variable s1, s2 Vectorised 2d cosine transform. Parameters ---------- f : function or method Function to apply 2D cosine trasnform to. f is called with f(x, y, *args). s1, s2 : 1d array Transformation variables. A grid of points will be made. args : tuple, optional Arguments to pass to f. m : int, optional Number of segments to break the integration interval into. Each segment will be between the zeros of the cos function, default m=20. ng : [2-20, 32, 64, 100], optional Number of gauss points to use in integration. Default ng=20. shanks_ind : int, optional Start position of intervals to start shanks extrapolation. default=None i.e. no extrapolation. Be careful when using shanks extrapolation; make sure you only begin to use it after the intgrand is well behaved. Returns ------- f : 1d array of float value of transform at s Notes ----- Careful with singularities. Because there is no way to increase the integration points at a particular sport the infinite behaviur may not be captured well. For example x**-0.5 should transform to sqrt(pi/2*w) but due to the sinularity at x=0 it does not converge well even using ng=100. """ #dims of array: # 0 or i dim is each of the 1st transform coordinate s1 # 1 or j dim is each integration inteval corresponding to s1, a1, b1 # 2 or k dim is each gauss point in interval a1, b1 # 3 or l dim is each of the 2nd transform coordinate s2 # 4 or m dim is each integration interval corresponding to s2, a2, b2 # 5 or n dim is each gauss point in interval a2, b2 # 5 dim will be summed to get integral of each interval a2, b2 # 4 dim will be summed to get shanks'ed to give transform at each s2 coord # 2 dim will be summed to get integral of each interval a1, b1 # 1 dim will be summed to get shanks'ed to give transform at each s1 coord si = np.atleast_1d(s1) sl = np.atleast_1d(s2) xk_, wk = gauss_legendre_abscissae_and_weights(ng) xn_, wn = xk_, wk # integration intervals zeros = np.zeros(m + 1, dtype=float) zeros[1:] = (2 * np.arange(m) + 1) * np.pi/2 aj = am = zeros[0:-1] bj = bm = zeros[1:] si = si[:, None, None, None, None, None] aj = aj[None, :, None, None, None, None] bj = bj[None, :, None, None, None, None] xk_ = xk_[None, None, :, None, None, None] wk = wk[None, None, :, None, None, None] aij = aj / si bij = bj / si bmaij = (bij - aij) / 2 # b minus a bpaij = (aij + bij) /2 # b plus a sl = sl[ None, None, None, :, None, None] am = am[None, None, None, None, :, None] bm = bm[None, None, None, None, :, None] xn_ = xn_[None, None, None, None, None, :] wn = wn[None, None, None, None, None, :] alm = am / sl blm = bm / sl bmalm = (blm - alm) / 2 # b minus a bpalm = (alm + blm) /2 # b plus a xijk = bmaij * xk_ + bpaij # xj_ are in [-1, 1] so need to transform to [a, b] xlmn = bmalm * xn_ + bpalm fijklmn=f(xijk, xlmn, *args) fijklmn*=np.cos(si*xijk) fijklmn*=np.cos(sl*xlmn) fijklmn*=wk fijklmn*=wn fijklmn*=bmaij#[:,:,0,:,:,:] fijklmn*=bmalm#[:,:,:,:,:,0] igral = np.sum(fijklmn, axis=5) if shanks_ind is None: igral = igral.sum(axis=4) else: #extrapolate igral.cumsum(axis=4 , out=igral) igral= np.apply_along_axis(shanks, 4, igral, shanks_ind) igral = np.sum(igral, axis=2) if shanks_ind is None: igral = igral.sum(axis=1) else: #extrapolate igral.cumsum(axis=1 , out=igral) igral= np.apply_along_axis(shanks, 1, igral, shanks_ind) return igral
##<Start of big commented block ## Below is a failed attempt at multiple transform. ## Instead see geotecha.mathematics.multi-transform #def cot(phi): # return 1/np.tan(phi) # #def csc(phi): # return 1.0/np.sin(phi) # #def vinv_laplace_2dcosine_transform(f, s1, s2, t, args=(), # m=20, ng=20, shanks_ind=None, # nlap=24, shift=0.0): # """Inverse laplace transform and 2D Cosine transform of f(x, y, t) # # Parameters # ---------- # f : function or method # function to apply cosine trasnform to. f is called with # f(x, y, *args) # s1, s2, t : 1d array # transformation variables. a grid of points will be made, s1 and s2 # are the x and y fourier directions, s3 is the t direction for the # inverse laplace # args : tuple, optional # arguments to pass to f # m : int, optional # number of segments to break the integration interval into. Each # segment will be between the zeros of the cos function, default=20 # ng : [2-20, 32, 64, 100], optional # number of gauss points to use in integration. # shanks_ind : int, optional # Start position of intervals to start shanks extrapolatoin. # default=None i.e. no extrapolation. # Be careful when using shanks extrapolation; make sure you only begin # to use it after the intgrand is well behaved. # nlap : even int, optional # number of integration points for talbot method inverse laplace. # if n is odd it will be rounded up to # nearest even number default n = 24 # shift : float # For inverse laplace transorm, shift contour to the right in case # there is a pole on the positive real axis. default shift=0.0 # # Returns # ------- # f : ndarray of shape (len(s1), len(s2), len(s3) # value of transform at s1, s2, s3 # # Notes # ----- # Careful with singularities in the fourier transform. Because there is # no way to increase the # integration points at a particular sport the infinite behaviur may not be # captured well. For example x**-0.5 should transform to sqrt(pi/2*w) but # due to the sinularity at x=0 it does not converge well even using ng=100. # # # # """ # # #dims of array: # # 0 or i dim is each of the 1st fourier transform coordinate s1 # # 1 or j dim is each integration inteval corresponding to s1, a1, b1 # # 2 or k dim is each gauss point in interval a1, b1 # # 3 or l dim is each of the 2nd fourier transform coordinate s2 # # 4 or m dim is each integration interval corresponding to s2, a2, b2 # # 5 or n dim is each gauss point in interval a2, b2 # # 6 or o dim is each of the inv laplace transform coordinate t # # 7 or p dim is the theta corresponding to t # # # 7 dim will be summed to get inv laplace for each t # # 5 dim will be summed to get integral of each interval a2, b2 # # 4 dim will be summed to get shanks'ed to give transform at each s2 coord # # 2 dim will be summed to get integral of each interval a1, b1 # # 1 dim will be summed to get shanks'ed to give transform at each s1 coord # # # to = np.atleast_1d(t) # nlap = nlap + nlap % 2 # if np.any(to==0): # raise ValueError('Inverse transform can not be calculated for t=0') # # Initiate the inv laplace stepsize # h = 2*np.pi/nlap; # theta = (-np.pi + (np.arange(nlap) + 1./2)*h) # # to = to[None, None, None, None, None, :, None] # thetap = theta[None, None, None, None, None, None, :] # zop = shift + nlap/to*(0.5017*thetap*cot(0.6407*thetap) - 0.6122 + 0.2645j*thetap) # dzop = nlap/to*(-0.5017*0.6407*thetap*(csc(0.6407*thetap)**2)+0.5017*cot(0.6407*thetap)+0.2645j) # ## theta = (-np.pi + (np.arange(self.n)+1./2)*h)[:, np.newaxis] ## z = self.shift + self.n/t*(0.5017*theta*cot(0.6407*theta) - 0.6122 + 0.2645j*theta) ## dz = self.n/t*(-0.5017*0.6407*theta*(csc(0.6407*theta)**2)+0.5017*cot(0.6407*theta)+0.2645j) ## inv_laplace = (np.exp(z * t) * self.f(z, *args) * dz).sum(axis=0) ## inv_laplace *= h / (2j * np.pi) # # # si = np.atleast_1d(s1) # sl = np.atleast_1d(s2) # # # # # xk_, wk = gauss_legendre_abscissae_and_weights(ng) # xn_, wn = xk_, wk # # # integration intervals # zeros = np.zeros(m + 1, dtype=float) # zeros[1:] = (2 * np.arange(m) + 1) * np.pi/2 # # aj = am = zeros[0:-1] # bj = bm = zeros[1:] # # si = si[:, None, None, None, None, None, None, None] # # aj = aj[None, :, None, None, None, None, None, None] # bj = bj[None, :, None, None, None, None, None, None] # # xk_ = xk_[None, None, :, None, None, None, None, None] # wk = wk[None, None, :, None, None, None, None, None] # # aij = aj / si # bij = bj / si # # bmaij = (bij - aij) / 2 # b minus a # bpaij = (aij + bij) /2 # b plus a # # # sl = sl[None, None, None, :, None, None, None, None ] # am = am[None, None, None, None, :, None, None, None] # bm = bm[None, None, None, None, :, None, None, None] # # xn_ = xn_[None, None, None, None, None, :, None, None] # wn = wn[None, None, None, None, None, :, None, None] # # alm = am / sl # blm = bm / sl # # bmalm = (blm - alm) / 2 # b minus a # bpalm = (alm + blm) /2 # b plus a # # # xijk = bmaij * xk_ + bpaij # xj_ are in [-1, 1] so need to transform to [a, b] # # xlmn = bmalm * xn_ + bpalm # # fijklmnop=f(xijk, xlmn, zop, *args) # fijklmnop*=np.cos(si*xijk) # fijklmnop*=np.cos(sl*xlmn) # # fijklmnop*=wk # fijklmnop*=wn # # fijklmnop*=bmaij#[:,:,0,:,:,:] # fijklmnop*=bmalm#[:,:,:,:,:,0] # # fijklmnop*=np.exp(zop * to) # fijklmnop*=dzop # fijklmnop*=h / (2j * np.pi) # # igral = np.sum(fijklmnop,axis=7) # igral = np.real(igral) # igral = np.sum(igral, axis=5) # ## igral = np.sum(fijklmn, axis=5) # # if shanks_ind is None: # igral = igral.sum(axis=4) # else: # #extrapolate # igral.cumsum(axis=4 , out=igral) # igral= np.apply_along_axis(shanks, 4, igral, shanks_ind) # # igral = np.sum(igral, axis=2) # # if shanks_ind is None: # igral = igral.sum(axis=1) # else: # #extrapolate # igral.cumsum(axis=1 , out=igral) # igral= np.apply_along_axis(shanks, 1, igral, shanks_ind) # # return igral # # # # # # # # # # # #def ilapcosine3(x, y, t, a, b, c): # """np.exp(-a * x) * np.exp(-b * y) / (1+t+c)""" # f= np.exp(-a * x) # f*=np.exp(-b * y) # f/=(1+t+c) # return f # #def ilapcosine3_(x, y, t, a, b, c): # return a / (a**2 + x**2) * b / (b**2 + y**2) *np.exp(-(c + 1) * t) # #class test_vinv_laplace_2dcosine_transform(unittest.TestCase): # """tests for vinv_laplace_2dcosine_transform""" # # def test_ilapcosine3(self): # s1 = np.array([0.5, 1, 1.6]) # s2 = np.array([0.6,1,2]) # t = np.array([0.2, 2, 4,]) # args=(1.2, 1.4, 0.8) # # shanks_ind=-5#None # nlap=24 # shift = 0.0 # assert_allclose(vinv_laplace_2dcosine_transform( # ilapcosine3, s1, s2, t, # args, shanks_ind=shanks_ind, # nlap=nlap, shift=shift), # ilapcosine3_(s1[:, None, None], # s2[None,:, None], # t[None, None, :], *args), atol=1e-8) ##>>END of big commented block. if __name__ == '__main__': import nose nose.runmodule(argv=['nose', '--verbosity=3', '--with-doctest']) # nose.runmodule(argv=['nose', '--verbosity=3'])