# Numerical integration in cython

I have a set of nested functions that I need to call multiple times. I know scipy.quad is pretty fast, but I will need to call the integrator recursively and want to remove as much overhead as possible.

I know I need to pre-allocate the arrays for some of the base functions at the bottom half of the code to get a speed up (the functions rvec, and below). However, my integration classes look like they can be optimized and I am puzzled as to how to go about this. Basically, I am preloading two vectors: weights and abscissas — which I will use for my Integrator class — these vectors are a constant and will never change. I have a class INorm that inherits from Integrator and a class klRatio that inherits from INorm. I think I need to use the __cinit__ function for the class initializations, but I am not sure the best way to go about this.

EDIT: I'm writing out Math Jax equations for all of the functions in a way that should match the code.

#!python
#cython: boundscheck=False
#cython: cdivision=True
#cython: wraparound=False

import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport sinh, cosh, sin, cos, sqrt, fabs, M_PI, floor, fmax, log

import h5py as h5

# C Math constants and functions
cdef double pi = float(3.141592653589793)

cdef extern from "complex.h":
double complex exp(double complex)

cdef extern from "complex.h":
double cabs(double complex)

cdef double[:] Abscissas2048
cdef double[:] Weights2048

with h5.File('LobattoNodes2048.h5','r') as hf:
Abscissas2048 = hf['Abscissas'][:]
Weights2048 = hf['Weights'][:]

# lobatto integrator
cdef class Integrator:
cdef double[:] abscissas
cdef double[:] weights
cdef int numEval
def __init__(self):
self.abscissas = Abscissas2048
self.weights = Weights2048
self.numEval = 2048
cdef double function(self,double x) except? 0.0:
raise NotImplementedError()
cdef double lobatto(self, double leftA, double rightB):
cdef double ax,dx,q
cdef int ii
ax = (leftA+rightB)/2.0
dx = (rightB-leftA)/2.0
cdef int startC = <unsigned int> self.numEval & <unsigned int> 0x1
cdef int endC = int((self.numEval+startC)/2)
q = startC*self.weights*self.function(ax)
for ii in range(startC,endC):
q += self.weights[ii]*(self.function(ax-dx*self.abscissas[ii])
+ self.function(ax+dx*self.abscissas[ii]))
return q*dx

cdef class INorm(Integrator):
cdef double leftB,rightB,Norm
cdef double[:,:] Pts
cdef double R,k,N
def __init__(self, double[:,:] Pts, double R, double k, int flag):
Integrator.__init__(self)
self.Pts = Pts
self.N = len(Pts)
self.R = R
self.k = k
self.leftB = 0.0
self.rightB = 2*pi
self.Norm = 1.0
if flag == 1:
self.updateNorm()
cdef double function(self, double theta) except? 0.0:
return self.R*Intensity(self.Pts, theta, self.R, self.k)/self.Norm
cdef double updateNorm(self):
self.Norm *= self.lobatto(self.leftB,self.rightB)
return 0.0

cdef class klRatio(INorm):
cdef double[:,:] Pts1,Pts2
cdef INorm N1
cdef INorm N2
cdef double KLoss, RingScore
def __init__(self, double[:,:] Pts1, double[:,:] Pts2, double R, double k):
INorm.__init__(self,Pts1, R, k,0)
self.N1 = INorm(Pts1,R,k,1)
self.N2 = INorm(Pts2,R,k,1)
self.KLoss = self.lobatto(self.leftB,self.rightB)
self.RingScore = fabs(log(self.N1.N/self.N2.N)) + self.KLoss
#print(self.RingScore)
cdef double function(self, double theta) except? 0.0:
return self.N1.function(theta)*log(self.N1.function(theta)
/self.N2.function(theta))
def get_RingScore(self):
return self.RingScore

# radial vectors of points from a position theta along the ring
cdef double[:] rvec(double[:,:] Pts, double theta, double R, Py_ssize_t N):
cdef int ii
cdef double [:] rvec = np.zeros(N)
for ii in range(0,N):
rvec[ii] = (R-Pts[ii,0]*cos(theta)
-Pts[ii,1]*sin(theta))
return rvec
# distance squared of points from a position theta along the ring
cdef double[:] r2Distance(double[:,:] Pts, double theta, double R, Py_ssize_t N):
cdef int ii
cdef double [:] r2 = np.zeros(N)
cdef double rx,ry
for ii in range(0,N):
rx = R*cos(theta)-Pts[ii,0]
ry = R*sin(theta)-Pts[ii,1]
r2[ii] = rx*rx + ry*ry
return r2

# Amplitude has to be normal vector to surface!
cdef double[:] Amplitude(double[:,:] Pts, double theta, double R, Py_ssize_t N):
cdef int ii
cdef double[:] numerator = rvec(Pts,theta,R,N)
cdef double[:] denominator = r2Distance(Pts,theta,R,N)
cdef double[:] Amp = np.zeros(N)
for ii in range(0,N):
Amp[ii] = sqrt(numerator[ii]/denominator[ii]/2/pi)
return Amp
# Phase returns an N vector of complex numbers
cdef double complex[:] Phase(double[:,:] Pts, double theta, double R, double k, Py_ssize_t N):
cdef int ii
cdef double[:] xComp = np.zeros(N)
cdef double[:] yComp = np.zeros(N)
cdef double complex[:] Pout = np.zeros(N,dtype=np.cdouble)
for ii in range(0,N):
xComp[ii] = (R*cos(theta)-Pts[ii,0])*cos(theta)
yComp[ii] = (R*sin(theta)-Pts[ii,1])*sin(theta)
Pout[ii] = exp(1j*k*(xComp[ii]+yComp[ii]))
return Pout
# Intensity, top level function to be integrated in classes
cdef double Intensity(double[:,:] Pts, double theta, double R, double k):
cdef size_t N
N = Pts.shape
cdef double[:] Amplitudes = Amplitude(Pts,theta,R,N)
cdef double complex[:] Phases = Phase(Pts,theta,R,k,N)
# pull out probability amplitudes as bras and kets
cdef double complex bra=0
cdef double complex ket
cdef int ii
cdef double Int
for ii in range(0,N):
bra += Amplitudes[ii]*Phases[ii]
ket = bra.conjugate()
Int = cabs(bra*ket)
return Int


The code is tested with something like this:

a = np.random.randn(30,2)/4
b = np.random.randn(30,2)/4
b[0:20,:] = a[0:20,:] # make some similarity among clusters

R = 1
k = 4*np.pi/R
zz = klRatio(a,b,R,k)
print(zz.get_RingScore)


OK, I'm going to stick with the scope of this code. Once I optimize this, I need to add another upper level so it's likely I'll be calling the quadrature operation a few hundred times, so getting it as fast as possible in Cython would be really nice. The theory specifically for the problem above:

The problem: There are two sets of 2-D coordinates, and we would like to measure how similar they are in a structural manner.

The solution: Solve for a modified KL-Divergence between two hypothetical radiation patterns given the provided coordinate information and measurement parameters.

We have two sets of coordinates $$\ P \$$ and $$\ Q \$$ with $$\ n \$$ and $$\ m \$$ member coordinates, respectively. There exists a normalizeable mapping function $$\ I(P,\theta,R,k) \$$ such that $$\ \int I(P,\theta,R,k) R d\theta \propto 1\$$. Here we have defined our measurement parameters in polar coordinates, where $$\ R \$$ is the radius of a ring shaped detector, $$\ \theta \$$ is a relative angle along the ring to some reference point, and k is the frequency of the radiation.

Since $$\ I() \$$ is normalizeable, we could effectively use the KL-divergence to a provide a quantitative divergence score. Since $$\ P \$$ and $$\ Q \$$ may have a different number of coordinates, the respective intensity measurements should scale accordingly to factor in this information. Without derivation, I include a symmetric penalty term for the loss function so that my modified KL-divergence, which I'll denote by the variable $$\ K \$$:

$$\ K = |\ln(n/m)| + \int R d\theta \tilde{I}(P,\theta,R,k) * \ln \frac{\tilde{I}(P,\theta,R,k)}{\tilde{I}(Q,\theta,R,k)} \$$

That is the final term that I wish to compute. To get to that term, I need to compute a normalized intensity: $$\ \tilde{I}(P,\theta,R,k)) = \frac{I(P,\theta,R,k)}{\int R d\theta I(P,\theta,R,k)}\$$. Since the intensity pattern is along a closed ring, the integral is closed at the end points ($$\0\$$ and $$\2 \pi \$$), I went with Lobatto integration. In the code above, I am starting with 2048 point integration but I will probably reduce the count when I incorporate the error estimate later.

The Intensity can be described from the direct measurement of probability amplitudes with sources defined by the $$\ P \$$ coordinates such that: $$\ I(P,\theta,R,k) = \Psi(P,\theta,R,k)^{\dagger} \Psi(P,\theta,R,k) \$$.

The probability amplitude given $$\P\$$ is defined as $$\ \Psi(P,\theta,R,k) = \sum^n_{i=1} A(P_i,\theta,R) \phi(P_i, \theta, R, k) \$$. The probability amplitude is complex, so it is described with an Amplitude $$\ A() \$$ and a phase $$\ \phi() \$$.

Where $$\ A(P_i,\theta,R) = \sqrt{\frac{R-x_i cos(\theta) - y_i sin(\theta)}{2 \pi ((R*cos(\theta)-x_i)^2 + (R*sin(\theta)-y_i)^2})} \$$. The amplitude term is the square root of the normal component $$\\hat{R}\$$ of the DC electric field in 2D. Here each coordinate member $$\ P_i = {x_i,y_i} \$$ is provided in the Cartesian system and is entirely enclosed within a radial distance $$\ R \$$. The code is not expected to work if the points are not enclosed.

The phase is defined as $$\ \phi(P_i,\theta,R,k) = \mathrm{exp}[\mathcal{i} k (R-x_i cos(\theta) - y_i sin(\theta)) ] \$$ The last part of the phase term is the dot product between the wave propagation vector and the normal component unit vector $$\ \hat{R} \$$ (normal to the circular ring).

That covers all the functions, more or less. My goal here is to be able to calculate several hundred kl-divergences in a reasonable amount of time as this is just one component in a larger algorithm. Thanks.

EDIT: In the future, I will have modified intensity/phase functions with an extra term and I will have to calculate a matrix of KL-divergences in order to deal with varying sub populations of structures. The number of integrations required will scales O(n^2) with this approach -- but dealing with 20 sub populations for each set, ~800 integrations is kind of a target I'd like to set for myself. The math is worked out for the sub population problem, I'm just trying to optimize all the simpler functions now before I re-work the larger problem. Thanks.

EDIT: For completeness, I'm adding the code for generating the lobatto weights and abscissas below:

import numpy as np
import h5py as h5

r8_epsilon = 2.220446049250313E-016

def r8vec_diff_norm_li ( n, a, b ):
value = 0.0
for i in range ( 0, n ):
value = max ( value, abs ( a[i] - b[i] ) )
return value

def r8vec_reverse ( n, a1 ):
a2 = np.zeros ( n )
for i in range ( 0, n ):
a2[i] = a1[n-1-i]
return a2

def getWeightsAbscissas(N):
A = np.zeros(N)
W = np.zeros(N)
if (N == 1):
A = -1.0
W = 2.0
return A, W
tol = 100 * r8_epsilon
for ii in range(0,N):
A[ii] = np.cos(np.pi * ii / (N-1))
Aold = np.zeros(N)
p = np.zeros([N,N])
while(True):
for ii in range(0,N):
Aold[ii] = A[ii]
for ii in range(0,N):
p[ii,0] = 1.0
p[ii,1] = A[ii]
for jj in range(2,N):
for ii in range(0,N):
p[ii,jj] = ((2*jj-1)*A[ii]*p[ii,jj-1]
+ (-jj+1)*p[ii,jj-2]) / jj
for ii in range(0,N):
A[ii] = Aold[ii] - (A[ii]*p[ii,N-1] - p[ii,N-2]) / (N*p[ii,N-1])
dif = r8vec_diff_norm_li(N,A,Aold)
if (dif <= tol): break
A = r8vec_reverse(N,A)
for ii in range(0,N):
W[ii] = 2.0/ (N*(N-1)*p[ii,N-1]**2)
return A,W

# write out weights and abscissas into an hdf5 file for pre-cached use
def WriteOutWA(N):
Ab,Wb = getWeightsAbscissas(N)
midPoint = int(np.floor(N/2))
Abscissas = Ab[midPoint:]
Weights = Wb[midPoint:]
# write out the h5 file
fileName = 'LobattoNodes{}.h5'.format(N)
h5f = h5.File(fileName,'w')
h5f.create_dataset("Abscissas", data=Abscissas)
h5f.create_dataset("Weights", data=Weights)
h5f.create_dataset("N", data=N)
h5f.close()

• Welcome to Code Review! You could probably make this a better question if you also add the describe the calculation that the code is performing with mathematical notation, using MathJax. – 200_success Jun 5 at 22:23