6
\$\begingroup\$

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[0]*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[0]
        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[0] = -1.0
        W[0] = 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()
\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – 200_success Jun 5 at 22:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.