def quadvgk(feval, fmin, fmax, tol1=<NUM_LIT>, tol2=<NUM_LIT>):
XK = np.array([-<NUM_LIT>, -<NUM_LIT>, -<NUM_LIT>, -<NUM_LIT>,<EOL>-<NUM_LIT>, -<NUM_LIT>, -<NUM_LIT>, <NUM_LIT:0.>,<EOL><NUM_LIT>, <NUM_LIT>, <NUM_LIT>,<EOL><NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>])<EOL>WK = np.array([<NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>,<EOL><NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>,<EOL><NUM_LIT>, <NUM_LIT>, <NUM_LIT>,<EOL><NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>])<EOL>WG = np.array([<NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>,<EOL><NUM_LIT>, <NUM_LIT>, <NUM_LIT>])<EOL>NK = WK.size<EOL>G = np.arange(<NUM_LIT:2>,NK,<NUM_LIT:2>)<EOL>tol1 = <NUM_LIT><EOL>tol2 = <NUM_LIT><EOL>Subs = np.array([[fmin],[fmax]])<EOL>NF = feval(np.zeros(<NUM_LIT:1>)).size<EOL>Q = np.zeros(NF)<EOL>neval = <NUM_LIT:0><EOL>while Subs.size > <NUM_LIT:0>:<EOL><INDENT>Subs = getSubs(Subs,XK)<EOL>M = (Subs[<NUM_LIT:1>,:] - Subs[<NUM_LIT:0>,:]) / <NUM_LIT:2><EOL>C = (Subs[<NUM_LIT:1>,:] + Subs[<NUM_LIT:0>,:]) / <NUM_LIT:2><EOL>NM = M.size<EOL>x = XK[:,None]*M + C<EOL>x = x.flatten()<EOL>FV = feval(x)<EOL>Q1 = np.zeros((NF, NM))<EOL>Q2 = np.zeros((NF, NM))<EOL>Q1 = np.dot(FV.reshape(NF, NK, NM).swapaxes(<NUM_LIT:2>,<NUM_LIT:1>),WK)*M<EOL>Q2 = np.dot(FV.reshape(NF, NK, NM).swapaxes(<NUM_LIT:2>,<NUM_LIT:1>)[:,:,<NUM_LIT:1>::<NUM_LIT:2>],WG)*M<EOL>ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)), <NUM_LIT:0>) < tol1 , (Subs[<NUM_LIT:1>,:] - Subs[<NUM_LIT:0>,:])  < tol2))[<NUM_LIT:0>]<EOL>Q = Q + np.sum(Q1[:,ind], axis=<NUM_LIT:1>)<EOL>Subs = np.delete(Subs, ind, axis=<NUM_LIT:1>)<EOL><DEDENT>return Q<EOL>
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: