4
\$\begingroup\$

I have a program doing a LOT of iteration (thousands to millions to hundreds of millions). It's starting to take quite a lot of time (few minutes, to a few days), and despite all my effort to optimize it, I'm still a bit stuck.

Profile: using cProfile via console call

ncalls     tottime  percall  cumtime  percall filename:lineno(function)
    500/1    0.018    0.000  119.860  119.860 {built-in method builtins.exec}
        1    0.006    0.006  119.860  119.860 Simulations_profiling.py:6(<module>)
      6/3    0.802    0.134  108.302   36.101 Simulations_profiling.py:702(optimization)
    38147    0.581    0.000  103.411    0.003 Simulations_profiling.py:270(overlap_duo_combination)
   107691   28.300    0.000  102.504    0.001 Simulations_profiling.py:225(duo_overlap)
 12615015   69.616    0.000   69.616    0.000 {built-in method builtins.round}

The first question, is what are the 2 first lines? I assumed that it was the program being called.

I'm going to replace the round() method by tolerance comparison in my if / else statements, thus avoiding this time consumption. I would like to optimize further the 2 following functions, but I can't find a new approach.

import itertools
import numpy as np

class Signal:
    def __init__(self, fq, t0, tf, width=0.3):
        self.fq = fq                                    # frequency in Hz
        self.width = float(width)                       # cathodic phase width in ms
        self.t0 = t0                                    # Instant of the first pulse in ms
        self.tf = tf                                    # End point in ms

        # List of instant at which a stim pulse is triggered in ms
        self.timeline = np.round(np.arange(t0, self.tf, 1/fq*1000), 3)

    def find_closest_t(self, t):
        val = min(self.timeline, key=lambda x:abs(x-t))
        id = np.where(self.timeline==val)[0][0]

        if val <= t or id == 0:
            return val, id
        else:
            return self.timeline[id-1], id-1

def duo_overlap(S1, S2, perc):

    pulse_id_t1, pulse_id_t2 = [], []

    start = max(S1.t0, S2.t0) - max(S1.width, S2.width)
    if start <= 0:
        start = 0
        start_id_S1 = 0
        start_id_S2 = 0
    else:
        start_id_S1 = S1.find_closest_t(start)[1]
        start_id_S2 = S2.find_closest_t(start)[1]

    stop = min(S1.tf, S2.tf) + max(S1.width, S2.width)

    for i in range(start_id_S1, len(S1.timeline)):
        if S1.timeline[i] > stop:
            break

        for j in range(start_id_S2, len(S2.timeline)):
            if S2.timeline[j] > stop:
                break

            d = round(abs(S2.timeline[j] - S1.timeline[i]), 3)  # Computation of the distance of the 2 point

            if S1.timeline[i] <= S2.timeline[j] and d < S1.width:
                pulse_id_t1.append(i)
                pulse_id_t2.append(j)
                continue

            elif S1.timeline[i] >= S2.timeline[j] and d < S2.width:
                pulse_id_t1.append(i)
                pulse_id_t2.append(j)
                continue

            else:
                continue

    return pulse_id_t1, pulse_id_t2

def overlap_duo_combination(signals, perc=0):

    overlap = dict()
    for i in range(len(signals)):
        overlap[i] = list()

    for subset in itertools.combinations(signals, 2):
        p1, p2 = duo_overlap(subset[0], subset[1], perc = perc)
        overlap[signals.index(subset[0])] += p1
        overlap[signals.index(subset[1])] += p2

    return overlap

Explanation of the program:

I have square signals of width Signal.width and of frequency Signal.fq starting at Signal.t0 and ending at Signal.tf. During the initialization of a Signal instance, the timeline is computed: it's a 1D-array of float in which each number corresponds to the instant at which a pulse is triggered.

Example:

IN: Signal(50, 0, 250).timeline
OUT: array([  0.,  20.,  40.,  60.,  80., 100., 120., 140., 160., 180., 200., 220., 240.])

A pulse is triggered at t = 0, t = 20, t = 40, ... Each pulse has a width of 0.3.

duo_overlap() takes 2 signals in input (and a percentage that we will keep fix at 0 for this example. This function computes the id of the the pulse for S1 and for S2 (ID in the timeline array) that overlap one with another.

Example:

If a pulse starts at t = 0 for S1 and a pulse starts at t = 0.2 for S2, since 0.2 - 0 = 0.2 < 0.3 (S1.width), the 2 pulses overlap.

I tried to optimize this function by looping only on the ID in which they can possibly overlap (start_id, stop) computed ahead.But as you can see, this function is still really time-consuming because of the high number of calls.

The last function, overlap_duo_combination() takes N signals in input as a list (or tuple / iterable) (2 <= N <= 6 in practice) and creates a dict() in which the key is the ID of the signal in the input list, and the value is a list of overlapping pulses ID (comparison 2 by 2 of the signals within the input list).

Example:

Input: signals = (S1, S2, S3) The pulse n°2 of S1 overlap with pulse n°3 of S2 and the pulse n°3 of S2 overlap with pulse n°5 of S3.

Output: dict[0] = [2] / dict[1] = [3, 3] / dict[2] = [5]

The 3 pops out twice for S2 because it will be add the first tile duo_overlap() is called on S1 and S2, and the second time when it is caleed on S2 and S3. I don't want to avoid duplicates since it's an information on how many different pulses are overlapping (in this case, 2 pulses are overlapping with the pulse n°3 of S2).

Would you have any idea, suggestion, code or whatever to redue the time complexity of this part of the code?

I am currently looking into PyCUDA implementation since I have a Nvidia 1080 Ti at disposal, but I don't know the C language. Would it be worth to switch to GPU this inner function that doesn't take long to execute when called but is called thousands of times?

Thanks for reading such a long post, and thanks for the help!

EDIT

Following the comment, numpy implementation of duo_overlap:

def duo_overlap_np(S1, S2, perc):
    p1_overlapping = np.zeros(shape = (len(S1.timeline)))
    p2_overlapping = np.zeros(shape = (len(S2.timeline)))

    start = max(S1.t0, S2.t0) - max(S1.width, S2.width)
    if start <= 0:
        start = 0
        start_id_S1 = 0
        start_id_S2 = 0
    else:
        start_id_S1 = S1.find_closest_t(start)[1]
        start_id_S2 = S2.find_closest_t(start)[1]

    stop = min(S1.tf, S2.tf) + max(S1.width, S2.width)

    for i in range(start_id_S1, len(S1.timeline)):
        if S1.timeline[i] > stop:
            break

        for j in range(start_id_S2, len(S2.timeline)):
            if S2.timeline[j] > stop:
                break

            d = round(abs(S2.timeline[j] - S1.timeline[i]), 3)  # Computation of the distance of the 2 point

            if S1.timeline[i] <= S2.timeline[j] and d < S1.width:
                p1_overlapping[i] = 1
                p2_overlapping[j] = 1
                continue

            elif S1.timeline[i] >= S2.timeline[j] and d < S2.width:
                p1_overlapping[i] = 1
                p2_overlapping[j] = 1
                continue

            else:
                continue

    return list(np.nonzero(p1_overlapping)[0]), list(np.nonzero(p2_overlapping)[0])

I'm quite sure I didn't get it right since it doesn't improve much. On 50 calls I get the following timing in second:

 Old: 8.20273494720459
 New: 8.188030004501343

EDIT2: In regards of the answer, the correct linspace() is:

import numpy as np

t0 = 0
tf = 300
fq = 60

t1 = np.round(np.arange(t0, tf, 1/fq*1000), 3)
t2 = np.round(np.linspace(t0, int(tf-(1/fq*1000)), int((tf-t0)*fq/1000)), 1)

And for the find_closest_t:

def find_closest_t(timeline, t):
    id = np.searchsorted(timeline, t)
    if id == 0:
        return timeline[0], id
    else:
        return timeline[id-1], id-1
\$\endgroup\$
  • \$\begingroup\$ I would consider avoiding the use of lists in duo_overlap and instead preallocate a numpy array with dimensions i,j using numpy.zeros. Elements that match your condition can then be set to 1 inside the loops. Afterwards, the indices containing nonzero elements found with numpy.nonzero() which should be equivalent to the lists you produce now (feel free to convert to a Python list at that point) \$\endgroup\$ – three_pineapples Mar 20 '18 at 10:13
  • \$\begingroup\$ @three_pineapples I just tried, see EDIT, but it didn't improve the performance :/ \$\endgroup\$ – Mathieu Mar 20 '18 at 10:52
  • \$\begingroup\$ In fact it even slowdowns the process compared to list. \$\endgroup\$ – Mathieu Mar 20 '18 at 13:48
  • \$\begingroup\$ Hmm, I guess maybe you aren't actually appending to the list that often? For a large number of appends I would have expected the numpy version to be faster. How big are those lists that get created? Also, it needed to be one 2D array which you set via myarray[i,j]=1 to reproduce the logic you had previously, not two 1D arrays (you overwrite elements in one of the arrays in each iteration of the inner loop). \$\endgroup\$ – three_pineapples Mar 20 '18 at 23:43
  • \$\begingroup\$ @three_pineapples The lists aren't that large. 10 to 20 items top. The issue is that the function is called 50K times or 5M times... Anyway, Niklas proposed a different approach using numpy that is really more efficient on SO: stackoverflow.com/questions/49380594/… \$\endgroup\$ – Mathieu Mar 21 '18 at 8:26
4
\$\begingroup\$

I only have a few minor things to offer

Having done some quick testing, np.linspace is significantly faster (about 60x) than np.arange

self.timeline = np.round(np.linspace(t0, self.tf, fq/1000), 3)

for find_closest_t, make use of the fact that self.timeline is in sorted order, since numpy has a function for that: https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.searchsorted.html This will make the search take O(log(n)) instead of O(n) where n is the size of timeline. You can use it to find the largest value less than or equal to t, and then just check the next value to find the smallest that is larger than t. The closest to t is one of the t

ind1 = np.searchsorted(self.timeline, t)
val1 = np.abs(self.timeline[ind1]-t)
ind2 = ind1+1
val2 = np.abs(self.timeline[ind2]-t)
ind, val = ind1, val1 if val1 < val2 else ind2, val2

You would need to add a conditional to check to make sure that ind1 wasn't the last index, or the one after it (in both situations, the last index should be val and ind)

if it weren't in sorted order, numpy's functions are optimized to be faster on arrays than those designed for any iterator:

val = np.min(np.abs(self.timeline-t))

does the exact same as your calculations for val, but faster

For your loops in duo_overlap_np, instead of doing S1.timeline[i] > stop, I would suggest:

s1_high = np.searchsorted(S1.timeline, stop, side='right')
s2_high = np.searchsorted(S2.timeline, stop, side='right')
for i in range(start_id_S1, s1_high):
    for j in range(start_id_S2, s2_high):

I don't think getting rid of the boolean check would save all that much time, though since for the S2 check you are running it over and over in the loop it could save a bit of time

lastly, on your return, converting numpy arrays to lists also takes a bit of computation time. Not sure what a simple way would be to change that up in this case though

\$\endgroup\$
  • \$\begingroup\$ Thanks for the complete reply. Although the changes are minor, even 1/100 of a second saved here reduce the overall computation of hours. I'll implement it this afternoon! \$\endgroup\$ – Mathieu Mar 20 '18 at 11:33
1
\$\begingroup\$

Code style

  • instead of iterating over the indices like this for j in range(start_id_S2, len(S2.timeline)): you can use enumerate like this: for i, s1 in enumerate(S1.timeline[start_id_S1:], start_id_S1):
  • there is no need for all the continue statements and the final else clause in your last if-part

performance

You can also use the fact that the pulses are in order to iterate over the signals concurrently. That way you wont have n! comparisons where you iterate over the whole tree

To do this, I made a small change to your Signal class:

class Signal:
    def __init__(self, freq, start, end, width=0.3):
        self.freq = freq                                    # frequency in Hz
        self.width = float(width)                       # cathodic phase width in ms
        self.start = start                                    # Instant of the first pulse in ms
        self.end = end                                    # End point in ms

        # List of instant at which a stim pulse is triggered in ms
        self.timeline = np.round(np.arange(start, end, 1000/freq), 3)
        self.pulses = np.stack((self.timeline, self.timeline + self.width), axis=1)

    def find_closest_t(self, t):
        val = min(self.timeline, key=lambda x:abs(x-t))
        id = np.where(self.timeline==val)[0][0]

        if val <= t or id == 0:
            return val, id
        else:
            return self.timeline[id-1], id-1

    def find_closest_t_np(self, t):
        idx = max(np.searchsorted(self.timeline, t) - 1, 0)
        return idx

    def __iter__(self):
        return iter(self.pulses)
        # or yield from map(tuple, self.pulses) # if you need tuples

    def __bool__(self):
        return bool(self.timeline.size)

To iterate over all the signals at the same time, we,

  1. assemble a dictionary with the signals and their first pulses.
  2. look for overlap in the pulses
  3. if this is a new overlap, yield it
  4. advance the iteration on the signal where the end of the pulse comes first
  5. if this iterator is exhausted, remove the signal from the dict with signals

code:

from collections import namedtuple

def duo_overlap_iter(signals, perc=0):
    pulse = namedtuple('Pulse', 'name iter index start end ')
    iters = ((i, iter(signal)) for i, signal in enumerate(signals) if signal)
    iters = {name: pulse(name, it, 0, *next(it)) for name, it in iters}
    seen = set()

    while iters:
        for overlap in find_overlap(iters.values()):
            if overlap not in seen:
                yield overlap
                seen.add(overlap)
        try:
            p0 = min(iters.values(), key=lambda x: (x.end, -x.start))
            iters[p0.name] = pulse(p0.name, p0.iter, p0.index + 1, *next(p0.iter))
        except StopIteration:
            del iters[p0.name]

To find the overlap, we use itertools.combinations. We yield an overlap as a frozenset with the name and index of the corresponding signals

def find_overlap(pulses):
    for p0, p1 in combinations(pulses, 2):
        p = frozenset(((p0.name, p0.index), (p1.name, p1.index)))
        if p1.start <= p0.end and p0.start <= p1.end:
            yield p

Sample Data

S0 = Signal(20 , 100, 0,)  # empty
S1 = Signal(50, 0, 250)
S2 = Signal(30, 10, 300, 2)
S3 = Signal(20, -10, 280, 2)
signals = S0, S1, S2, S3

sample result

list(duo_overlap_iter(signals))
[frozenset({(1, 2), (3, 1)}), frozenset({(3, 3), (1, 7)}), frozenset({(1, 12), (3, 5)})]

Final results

To get the results in the way your code presents it, you can do something like this:

def overlap_duo_comb_iter(signals, perc=0):

    overlap = {i: [] for i, _ in enumerate(signals)}

    for (s0, i0), (s1, i1) in duo_overlap_iter(signals):
        overlap[s0].append(i0)
        overlap[s1].append(i1)

    return overlap

numpy code

During my revision of the timing, I noticed in your edited code using numpy, you still don't use the fact the signals are sorted. You iterate over the whole timeline instead of stopping once the 2nd signal can not overlap anymore.

def duo_overlap_np(S1, S2, perc):
    p1_overlapping = np.zeros_like(S1.timeline)
    p2_overlapping = np.zeros_like(S2.timeline)

    start = max(S1.start, S2.start)
    start_id_S1 = S1.find_closest_t_np(start)

    stop = min(S1.pulses[-1][1], S2.pulses[-1][1])
    for i, (s1, s1_end) in enumerate(S1.pulses[start_id_S1:], start_id_S1):
        if s1 > stop:
            break

        start_id_S2 = S2.find_closest_t_np(s1)
        for j, (s2, s2_end) in enumerate(S2.pulses[start_id_S2:], start_id_S2):
            if s2 > s1_end:
                break
            if s1 > s2_end:
                continue
            p1_overlapping[i] = 1
            p2_overlapping[j] = 1

    return list(np.nonzero(p1_overlapping)[0]), list(np.nonzero(p2_overlapping)[0])

Timings

print(overlap_duo_combination(signals))
% timeit overlap_duo_combination(signals)
{0: [], 1: [2, 7, 12], 2: [], 3: [1, 3, 5]}
1.33 ms ± 72.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
print(overlap_duo_combination(signals, func=duo_overlap_np))
assert overlap_duo_combination(signals) == overlap_duo_combination(signals, func=duo_overlap_np)
% timeit overlap_duo_combination(signals, func=duo_overlap_np)
{0: [], 1: [2, 7, 12], 2: [], 3: [1, 3, 5]}
267 µs ± 4.75 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
print(overlap_duo_comb_iter(signals))
assert overlap_duo_combination(signals) == overlap_duo_comb_iter(signals,)
% timeit overlap_duo_comb_iter(signals)
{0: [], 1: [2, 7, 12], 2: [], 3: [1, 3, 5]}
600 µs ± 12.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
print(list(duo_overlap_iter(signals)))
% timeit list(duo_overlap_iter(signals))
[frozenset({(1, 2), (3, 1)}), frozenset({(3, 3), (1, 7)}), frozenset({(1, 12), (3, 5)})]
605 µs ± 33.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Code can all be found on github

Conclusion

The adapted numpy code is the fasted for this set of signals, but only works with fixed-length pulses. So if your pulses are fixed length, use the numpy code, if you just get pulses, with a start and end, my code can handle those too

\$\endgroup\$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.