# Optimize Cython code with np.ndarray contained

Can someone help me further optimize the following Cython code snippets? Specifically, a and b are np.ndarray with int value (range(256)) in them. They are one dimension arrays with dynamic length. resultHamming is a one-dimension array with float value in it (dynamic length). bits is an int list (size 256).

The function is to compare two dynamic length bit vector, and return a similarity value as the distance of the two, where the length of each vector is a multiple of 2048-bit (256 bytes). I want to find the best match between these two bit vector by comparing each 2048-bit block, where each bit vector is represented as ndarray (read the bit sequence byte by byte, thus each position is range from 0 to 2^8 = 256). Rule for matching is to find global minimum distance between all block pairs and allow one block in A to be matched with more than one block in B if they have smaller distance. Always compare the smaller size vector against the larger one.

The following code assumes b vector is smaller. We can limit resultHamming to be smaller than size of numArrayB and only record numArrayB smallest distance value, but need to track the current size when inserting new value into it. Even with current case (record all the pairwise distance), we actually know the final size of resultHamming at the beginning.

def compare(a, b):
cdef double dis, hammingTotal = 0
cdef int numArrayA = int(a.size/256)
cdef int numArrayB = int(b.size/256)
cdef int i, j, k, l, index
bits = list(xrange(256))
# Prepare a bit number table for fast query
for l in xrange(256):
# nnz() counts the number of 1s in value
bits[l] = nnz(l)

resultHamming = []
for i in xrange(numArrayB):
# Count the number of 1-bits in i-th block of B
onesB = sum(bits[b[k+256*i]] for k in xrange(256))
for j in xrange(numArrayA):
# Count the number of 1-bits in j-th block of A
onesA = sum(bits[a[k+256*j]] for k in xrange(256))
# Calculate the hamming distance between i-th block of B and j-th block of A
hammingCur = sum(bits[b[i*256+k] ^ a[j*256+k]] for k in xrange(256))
dis = (hammingCur) / (onesA + onesB)
# Insertion current dis to resultHamming with sorted order
k = len(resultHamming) - 1
if dis >= resultHamming[-1]:
resultHamming.append(dis)
else:
resulthamming.append(resultHamming[k])
while k > 0 and resultHamming[k-1] > dis:
resultHamming[k] = resultHamming[k-1]
k -= 1
resultHamming[k] = dis

# Extract k smallest distance from the distance array
for k in xrange(numArrayB):
hammingTotal += resultHamming[k]
return round(hammingTotal/(numArrayB), 3)

• What are you trying to achieve here? What is the specification of your function compare? – Gareth Rees Mar 6 '14 at 20:13
• Update the specification and code a little bit. Previously, the rule does not allow one block in A to be matched with several blocks in B, thus it would be even more tricky to find the final k smallest distance value. – Rain Lee Mar 6 '14 at 22:13
• The way you have put it, distance from A to B differs from the distance from B to A. What are you looking for? In addition, I wonder, why are you normalizing hamming distance with the number of set bits in each block? As a general rule for optimizing with cython, try to avoid python function calls if possible, i.e. replace pythonic sum with numpy sum, replace bit number table with numpy array, etc... – Blaz Bratanic Mar 6 '14 at 23:11
• Compare(A, B) tries to find the best block by block match between A and B. The above code assume B.size is smaller than A.size, there is an extra step (I omit for space) to switch A and B if B.size > A.size, or add another block to put A in outer layer for loop (Always put the smaller size input in outer layer for loop, and compare the smaller one against the larger one). Since I am choosing dis values globally from the pairwise distance matrix, compare(A,B) should equal to compare(B,A). The flaw exists in answer 1 though. The intention and problem of normalizing is covered in comments below. – Rain Lee Mar 7 '14 at 6:30

As a general rule for optimizing with cython, try to avoid python function calls if possible, i.e. replace pythonic sum with numpy sum or explicit for loop, replace bit number table with numpy array, and try to declare types for all variables.

You can also precalculate most of the variables used in distance calculation and use a separate array to hold minimum distances for each block. I quickly put something together and got the code below.

Just wondering, what is the rationale for normalizing the distance with the number of set bits?

import cython
import numpy as np
cimport numpy as np

cdef bitCount(unsigned char value):
cdef count = 0
while(value):
value &= value - 1
count += 1
return(count)

def compare(np.ndarray[np.uint8_t, ndim=1, mode='c'] a,
np.ndarray[np.uint8_t, ndim=1, mode='c'] b):

cdef int nblocksA = int(a.size/256)
cdef int nblocksB = int(b.size/256)
cdef int i, j, k

cdef double partial_hamming_distance

# Prepare a bit number table for fast query
cdef np.ndarray[np.uint8_t, ndim=1] bits = np.zeros(256, dtype=np.uint8)
for i in range(256):
bits[i] = bitCount(i)

cdef np.ndarray[np.float_t, ndim=1] min_distances = np.ones(nblocksB,
dtype=np.float)
min_distances *= 1000000

cdef np.ndarray[np.float_t, ndim=1] set_bitsA = np.zeros(nblocksA, dtype=np.float)
for i in range(nblocksA):
for k in range(256):
set_bitsA[i] += bits[a[k + 256 * i]]

cdef np.ndarray[np.float_t, ndim=1] set_bitsB = np.zeros(nblocksB, dtype=np.float)
for i in range(nblocksB):
for k in range(256):
set_bitsB[i] += bits[b[k + 256 * i]]

for i in range(nblocksB):
for j in range(nblocksA):
partial_hamming_distance = 0

for k in range(256):
partial_hamming_distance += (bits[b[i*256+k] ^ a[j*256+k]]
/ (set_bitsB[i] + set_bitsA[j]))

if partial_hamming_distance < min_distances[i]:
min_distances[i] = partial_hamming_distance

return np.sum(min_distances)

• Thanks for the beautiful code style and cython advice. You definitely got my intention. 1) Does it better if use np.empty for ndarray initialization? 2) You have three 2+ layer for loops, what about merging them into one? 3) Your partial hamming_distance computation logic is my original version. There might be a flaw: assume A=(a1,a2), B=(b1,b2). And the pairwise distance matrix is: (b1,a1)=46, (b1,a2)=64, (b2,a1)=78, (b2,a2)=85. The above code would return (46+87)/2=62, b1,b2 both match with a1; but a better match would be return (46+64)/2=55, a1,a2 both match b1, compare(a,b)!=compare(b,a). – Rain Lee Mar 7 '14 at 5:52
• Sorry I failed to put all my words in above comments. The distance got from above is not global minimum. And The reason for normalizing the distance like is mainly to return the final distance in scale of 1 instead of 2048. The above code is actually from a fuzzy hashing algorithm called mvHash-B. It might be problematic in case where (set_bitsB[i]+set_bitsA[j]) > 2048 since each block would have 2048 bits at most. BTW, the parameter type of function bitCount() is unsigned char, what's the range of "value" here? – Rain Lee Mar 7 '14 at 7:05
• @RainLee 1. np.empty is even better choice and will be marginally faster than initializing values with zero. 2. the first two for loops precalculate number of set bits in each block. And the third one does the actual calculation. You might merge second for loop with the actual calculation since its iterating over the same range. 3. I will look into it when I have some more time. – Blaz Bratanic Mar 7 '14 at 8:46
• @RainLee Oh, unsigned char represents an unsigned 8bit value [0, 256) – Blaz Bratanic Mar 7 '14 at 8:53
• @RainLee Also... you can compile cython with cython -a MODULE.py. This will generate html showing each line of source code colored white through various shades of yellow. The darker the yellow color, the more dynamic Python behaviour is still being performed on that line. For each line that contains some yellow, you need to add more static typing declarations. – Blaz Bratanic Mar 7 '14 at 13:58