Source code for GPy.util.quad_integrate

"""
The file for utilities related to integration by quadrature methods
- will contain implementation for gaussian-kronrod integration.

"""
import numpy as np

[docs]def getSubs(Subs, XK, NK=1): M = (Subs[1, :] - Subs[0, :]) / 2 C = (Subs[1, :] + Subs[0, :]) / 2 I = XK[:, None] * M + np.ones((NK, 1)) * C # A = [Subs(1,:); I] A = np.vstack((Subs[0, :], I)) # B = [I;Subs(2,:)] B = np.vstack((I, Subs[1, :])) # Subs = [reshape(A, 1, []); A = A.flatten() # reshape(B, 1, [])]; B = B.flatten() Subs = np.vstack((A,B)) # Subs = np.concatenate((A, B), axis=0) return Subs
[docs]def quadvgk(feval, fmin, fmax, tol1=1e-5, tol2=1e-5): """ numpy implementation makes use of the code here: http://se.mathworks.com/matlabcentral/fileexchange/18801-quadvgk We here use gaussian kronrod integration already used in gpstuff for evaluating one dimensional integrals. This is vectorised quadrature which means that several functions can be evaluated at the same time over a grid of points. :param f: :param fmin: :param fmax: :param difftol: :return: """ XK = np.array([-0.991455371120813, -0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691, -0.405845151377397, -0.207784955007898, 0., 0.207784955007898, 0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]) WK = np.array([0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529]) # 7-point Gaussian weightings WG = np.array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277, 0.129484966168870]) NK = WK.size G = np.arange(2,NK,2) tol1 = 1e-4 tol2 = 1e-4 Subs = np.array([[fmin],[fmax]]) # number of functions to evaluate in the feval vector of functions. NF = feval(np.zeros(1)).size Q = np.zeros(NF) neval = 0 while Subs.size > 0: Subs = getSubs(Subs,XK) M = (Subs[1,:] - Subs[0,:]) / 2 C = (Subs[1,:] + Subs[0,:]) / 2 # NM = length(M); NM = M.size # x = reshape(XK * M + ones(NK, 1) * C, 1, []); x = XK[:,None]*M + C x = x.flatten() FV = feval(x) # FV = FV[:,None] Q1 = np.zeros((NF, NM)) Q2 = np.zeros((NF, NM)) # for n=1:NF # F = reshape(FV(n,:), NK, []); # Q1(n,:) = M. * sum((WK * ones(1, NM)). * F); # Q2(n,:) = M. * sum((WG * ones(1, NM)). * F(G,:)); # end # for i in range(NF): # F = FV # F = F.reshape((NK,-1)) # temp_mat = np.sum(np.multiply(WK[:,None]*np.ones((1,NM)), F),axis=0) # Q1[i,:] = np.multiply(M, temp_mat) # temp_mat = np.sum(np.multiply(WG[:,None]*np.ones((1, NM)), F[G-1,:]), axis=0) # Q2[i,:] = np.multiply(M, temp_mat) # ind = np.where(np.logical_or(np.max(np.abs(Q1 -Q2) / Q1) < tol1, (Subs[1,:] - Subs[0,:]) <= tol2) > 0)[0] # Q = Q + np.sum(Q1[:,ind], axis=1) # np.delete(Subs, ind,axis=1) Q1 = np.dot(FV.reshape(NF, NK, NM).swapaxes(2,1),WK)*M Q2 = np.dot(FV.reshape(NF, NK, NM).swapaxes(2,1)[:,:,1::2],WG)*M #ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)/Q1), 0) < difftol , M < xtol))[0] ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)), 0) < tol1 , (Subs[1,:] - Subs[0,:]) < tol2))[0] Q = Q + np.sum(Q1[:,ind], axis=1) Subs = np.delete(Subs, ind, axis=1) return Q
[docs]def quadgk_int(f, fmin=-np.inf, fmax=np.inf, difftol=0.1): """ Integrate f from fmin to fmax, do integration by substitution x = r / (1-r**2) when r goes from -1 to 1 , x goes from -inf to inf. the interval for quadgk function is from -1 to +1, so we transform the space from (-inf,inf) to (-1,1) :param f: :param fmin: :param fmax: :param difftol: :return: """ difftol = 1e-4 def trans_func(r): r2 = np.square(r) x = r / (1-r2) dx_dr = (1 + r2)/(1-r2)**2 return f(x)*dx_dr integrand = quadvgk(trans_func, -1., 1., difftol, difftol) return integrand