# Optimising compatibility-conflict graph solution - python recursive routine

### Background

I'm implementing an algorithm which localises overlapping sound sources using the time-difference-of-arrivals across 3 sensors [1]. Compatible 'triples' (e.g. with only 2 common sensors) need to be joined up to make larger N sensor graphs. The compatible triples are found using a recursive routine (combine_all) until no new solutions are found.

The routine runs well when there are small # of triples, but I've empirically observed for every 10X increase in triple number, there's a ~10^5.5 increase in run-time (using linear regression) - which is problematic (meaning up to 12 hours of run-time on my actual data).

##### Reference

[1] Kreissig & Yang 2013, Fast and reliable TDOA assignment in multi-source reverberant environments, ICASSP 2013 paper link.

### Code and initial profiling

The combine_all routine accepts the compatibility-conflict matrix describing the compatibility/conflict of all triple pairs. From the currently available triples - their in/compatibility to the current solution are checked and kept or eliminated. The compatible nodes are checked by get_Nvl, and incompatible nodes are checked by get_not_Nvl.

Initial profiling with iPython %lprun told me the get_Nvl and get_not_Nvl is where ~80% (40% and 40% each) is spent in combine_all - and I've tried my best to optimise the current code with no luck.

import numpy as np
from itertools import chain, product

def combine_all(Acc, V, l, X):
'''
Parameters
----------
Acc : (N_triples,N_triples) np.array
The compatibility-conflict graph. Value of 1 means compatible node pair, -1 is incompatible, 0 is undefined.
V_t : set
V_tilda. The currently considered vertices (a sub-set of V, all vertices)
l : set
The solution consisting of the currently compatible vertices.

Returns
-------
solutions_l : list with sublists
A somewhat messy output - must be run through format_combineall to get
nice output.

'''
# determine N_v(l) and !N_v(l)
# !N_v(l) are the vertices incompatible with the current solution
N_vl = get_Nvl(Acc, V, l)
N_not_vl = get_NOT_Nvl(Acc, V, l)

solutions_l = []
if len(N_vl) == 0:
solutions_l.append(l)

else:
# remove conflicting neighbour
V = V.difference(N_not_vl)
# unvisited compatible neighbours
Nvl_wo_X = N_vl.difference(X)
for n in Nvl_wo_X:
Vx = V.difference(set([n]))
lx = l.union(set([n]))
solution = combine_all(Acc, Vx, lx, X)
if solution:
solutions_l.append(solution)
X = X.union(set([n]))
return solutions_l

def get_Nvl(Acc, V_t, l):
'''Checks for compatible vertices

Essentially, two for loops run across the
available triples (v in V_t) and   the current solution set (u in l)
- to access the Acc[v,u] entries. If all Acc[v,u] values
are +1, then v is compatible with the current solution l.
If any of the Acc[v,u] values is -1  then that v is not compatible anymore.

Returns
-------
Nvl : set
Solution of vertices that are compatible to at least one other vertex
and not in conflict with any of the other vertices.
'''
Nvl = []
if len(l)>0:
for v in V_t:
for u in l:
if Acc[v,u]==1:
Nvl.append(v)
elif Acc[v,u]==-1:
if v in Nvl:
Nvl.pop(Nvl.index(v))
return set(Nvl)
else:
return V_t

def get_Nvl_fast(Acc, V_t, l):
'''See CombineAll for docs
'''
if len(l)>0:
all_uv = np.array(np.meshgrid(V_t, l)).T.reshape(-1,2)
def get_acc(X):
return Acc[X[0], X[1]]
Acc_values = np.apply_along_axis(get_acc, 1, all_uv)
rows_w_min1 = np.where(Acc_values<0)
v_vals_w_conflicts = np.unique(all_uv[rows_w_min1,0])
Nvl = np.setdiff1d(V_t, v_vals_w_conflicts)
return Nvl
else:
return V_t

def get_NOT_Nvl(Acc:np.array, V:set, l:set):
N_not_vl = []
if len(l)>0:
for v in V:
for u in l:
if Acc[v,u]==-1:
N_not_vl.append(v)
elif Acc[v,u]==1:
if v in N_not_vl:
N_not_vl.pop(N_not_vl.index(v))
else:
N_not_vl = []
return set(N_not_vl)

# ---- not performance related section -- formats output into easily readable
# form
def format_combineall(output):
semiflat = flatten_combine_all(output)
only_sets = []
for each in semiflat:
if isinstance(each, list):
for every in each:
if isinstance(every, set):
only_sets.append(every)
elif isinstance(each, set):
only_sets.append(each)
return only_sets

def flatten_combine_all(entry):
if isinstance(entry, list):
if len(entry)==1:
return flatten_combine_all(entry[0])
else:
return list(map(flatten_combine_all, entry))
elif isinstance(entry, set):
return entry
else:
raise ValueError(f'{entry} can only be set or list')

if __name__ == '__main__':
# compatibility-conflict graph from [1]
A = np.array([[ 0, 1, 0, 0,-1,-1],
[ 1, 0, 1, 1, 0, 1],
[ 0, 1, 0,-1, 1, 0],
[ 0, 1,-1, 0,-1, 0],
[-1, 0, 1,-1, 0, 1],
[-1, 1, 0, 0, 1, 0]])

qq = combine_all(A, set(range(6)), set([]), set([]))
neat_output = format_combineall(qq)
# >>> print(neat_output)
# >>> [{0, 1, 2}, {0, 1, 3}, {1, 2, 4, 5}, {1, 3, 5}]


### No luck with optimisation experiments

Since the get_Nvl + get_not_Nvl is where most of the time spent - I've focussed on performing optimisations there. To begin with I worked only get_Nvl (both Nvl functions have very similar structure. In the current get_Nvl implementation - there are two for loops, with i,j referencing values in the Acc compatibility-conflict matrix. I've tried multiple things that have lead to no improvement or even increase in runtime, of which here I report the two that I can now recollect.

1. Converting serial loop-in-loop (i,j) into a direct i,j product - that can be used to check the values of Acc in a map call.
2. converting the loop-in-loop into a numpy iterable form (np.apply_along_axis) (e.g. see get_Nvl_fast)
3. Converting the if,else flow into dictionary calls. (e.g. action_to_take[Acc[v,u]] instead of if Acc[v,u]==-1:...)

### Solutions to speeding up the code?

I've run out of ideas on how to speed up the code in Python - and am now considering using a C++ implementation that I can call from Python (though my C++ is non-existent/rusty).

I'd be grateful for ideas on how to speed up my Python code before investing too deep into another language!

## Style & Maintenance

Prior to optimisation fussiness,

Use PEP484 type hints.

Docstrings are typically in double quotes rather than single quotes.

Do a PEP8 pass with a linter or good IDE; they will point out e.g. that you need spaces around ==.

You're so close to having a unit test! Write an assert.

## Performance

Some of this will make a difference and some won't. It's good that you're profiling; keep doing that.

You've missed some opportunities for early return; pursue this, for example when you get_NOT_Nvl(Acc, V, l).

Rather than V = V.difference, just V -=.

Rather than set([n]), prefer {n}; but better yet don't do a set-to-set operation at all. Clone the set and then use .discard and .add as needed, since this is a single-element operation.

get_Nvl needs to operate on Nvl as a set instead of as a list cast to a set. I anticipate this helping, since Nvl.pop() and Nvl.index() are much more inefficient than a set discard. The same applies to get_NOT_Nvl.

More broadly, Python is bad at recursion. If you absolutely must have a recursive solution, as you fear: it is time to brush up on C++.

## Suggested

import numpy as np

def combine_all(Acc: np.ndarray, V: set[int], l: set[int], X: set[int]) -> list[set[int]]:
"""
Parameters
----------
Acc : (N_triples,N_triples) np.array
The compatibility-conflict graph. Value of 1 means compatible node pair, -1 is incompatible, 0 is undefined.
V_t : set
V_tilda. The currently considered vertices (a sub-set of V, all vertices)
l : set
The solution consisting of the currently compatible vertices.

Returns
-------
solutions_l : list with sublists
A somewhat messy output - must be run through format_combineall to get
nice output.

"""
# determine N_v(l) and !N_v(l)
# !N_v(l) are the vertices incompatible with the current solution
N_vl = get_Nvl(Acc, V, l)
if len(N_vl) == 0:
return [l]

solutions_l = []
# remove conflicting neighbour
V -= get_NOT_Nvl(Acc, V, l)
# unvisited compatible neighbours
N_vl -= X
X = set(X)

for n in N_vl:
Vx = set(V)
lx = set(l)
solution = combine_all(Acc, Vx, lx, X)
if solution:
solutions_l.append(solution)
return solutions_l

def get_Nvl(Acc: np.ndarray, V_t: set[int], l: set[int]) -> set[int]:
"""Checks for compatible vertices

Essentially, two for loops run across the
available triples (v in V_t) and   the current solution set (u in l)
- to access the Acc[v,u] entries. If all Acc[v,u] values
are +1, then v is compatible with the current solution l.
If any of the Acc[v,u] values is -1  then that v is not compatible anymore.

Returns
-------
Nvl : set
Solution of vertices that are compatible to at least one other vertex
and not in conflict with any of the other vertices.
"""
if len(l) < 1:
return V_t

Nvl = set()
for v in V_t:
for u in l:
a = Acc[v, u]
if a == 1:
elif a == -1:
return Nvl

def get_NOT_Nvl(Acc: np.array, V: set[int], l: set[int]) -> set[int]:
if len(l) < 1:
return set()

N_not_vl = set()
for v in V:
for u in l:
if Acc[v, u] == -1:
elif Acc[v, u] == 1:
return N_not_vl

def format_combineall(output: list) -> list[set[int]]:
"""formats output into easily readable
form"""
semiflat = flatten_combine_all(output)
only_sets = []
for each in semiflat:
if isinstance(each, list):
for every in each:
if isinstance(every, set):
only_sets.append(every)
elif isinstance(each, set):
only_sets.append(each)
return only_sets

def flatten_combine_all(entry: list) -> list:
if isinstance(entry, list):
if len(entry) == 1:
return flatten_combine_all(entry[0])
else:
return list(map(flatten_combine_all, entry))
if isinstance(entry, set):
return entry
raise TypeError(f'{entry} can only be set or list')

def test() -> None:
# compatibility-conflict graph from [1]
A = np.array([[ 0, 1, 0, 0,-1,-1],
[ 1, 0, 1, 1, 0, 1],
[ 0, 1, 0,-1, 1, 0],
[ 0, 1,-1, 0,-1, 0],
[-1, 0, 1,-1, 0, 1],
[-1, 1, 0, 0, 1, 0]])

qq = combine_all(A, set(range(6)), set(), set())
neat_output = format_combineall(qq)
assert neat_output == [{0, 1, 2}, {0, 1, 3}, {1, 2, 4, 5}, {1, 3, 5}]

if __name__ == '__main__':
test()
`
• Thanks a lot for the suggestions! Set based formulation is much more elegant - however, as you predicted, sadly not much faster :|. Need to go down the Cpp path I guess.. Aug 9, 2022 at 6:25