Timeout functionality

[Jnelen]

Hi there

First of all I'd like to say that I think this is a really nice and useful package, so thanks for developing it!

However, I encountered some performance bottleneck when calculating the symmetry-corrected RMSD for larger molecules. I was performing docking calculations using Drugbank molecules, and these molecule in particular gave problems:
https://go.drugbank.com/drugs/DB14981
https://go.drugbank.com/drugs/DB01666
https://go.drugbank.com/drugs/DB05808
I am assuming this is due to large molecule size, and a potential combinatorial explosion when trying to figure out the optimal RMSD when considering the symmetry options.
Maybe it could be interesting to add some kind of time-out functionality for these cases? For standard settings, there probably should be no timeout to keep default behavior, but it could be a useful feature. For example if you set the timeout as 10, it will raise an error or return NaN if it doesn't return anything after running for 10 seconds. I could try to take a stab at implementing something like this and make a pull request, unless you don't think this is a valuable feature?

[RMeli]

Hello @Jnelen, thank you for your report.

I am assuming this is due to large molecule size, and a potential combinatorial explosion when trying to figure out the optimal RMSD when considering the symmetry options.

Looking at the molecules, I think the problem is the very high number of graph isomorphisms rather than the size.

I have encountered similar issues with a few similar molecules to the ones you reported. Out of curiosity, which backend are you using? I noticed that networkx had troubles with such molecules, while graph-tool remained reasonably fast. If you used the former, would you mind trying to install graph-tool and see if it can handle them? (If graph-tool is installed, it will supersede networkx; I still need to address #68...)

Maybe it could be interesting to add some kind of time-out functionality for these cases?

I think this might be useful, especially if

• The default is to have no timeout (as you suggest)
• The timeout raises an error rather than returning a NaN
  The latter is my personal opinion, so it's up for discussion if you have good arguments about returning NaNs.

[RMeli]

As for backend, is there a way to check which one is being used?

Unfortunately, not really. That would be a good improvement, but the real solution should be #68. Unfortunately I had no time to implement this yet.

I'm running off a singularity image and installed graph-tool using conda install -c conda-forge graph-tool, but it would be good to make sure it's properly recognized by spyrmsd.

If you have it installed, it will be automatically used before networkx.

---

I just had a quick try with the first molecule (on macOS M1), and I get

$ time python -m spyrmsd DB14981.sdf DB14981.sdf
0.00000
python -m spyrmsd DB14981.sdf DB14981.sdf  99.75s user 0.34s system 94% cpu 1:46.46 total

with networkx and

$ time python -m spyrmsd DB14981.sdf DB14981.sdf
time python -m spyrmsd test.sdf test.sdf
0.00000
python -m spyrmsd DB14981.sdf DB14981.sdf  7.77s user 0.28s system 99% cpu 8.096 total

with graph-tool. Quite long, but still reasonable.

However, for the last molecule I have

$ time python -m spyrmsd DB05808.sdf DB05808.sdf
[2]    41710 killed     python -m spyrmsd test2.sdf test2.sdf
python -m spyrmsd DB05808.sdf DB05808.sdf  223.30s user 36.83s system 87% cpu 4:57.13 total

even with graph-tool. So I agree that a timeout would be helpful. But maybe a simpler alternative would just be to use the timeout command (on Linux/macOS; I don't know much about Windows)?

timeout 20 python -m spyrmsd test.sdf test.sdf

(For macOS, one needs to install coreutils with brew install coreutils and then use gtimeout.)

I'm also wondering if there is some graph metric that can single out these pathological cases (eventually to automatically turn on a timeout). I'll look into it.

[Jnelen]

Hi @RMeli,

I'm using spyrmsd to calculate many RMSD values for different molecules as part of a for loop. As you indicated, the main problem was the last molecule, actually resulting in the whole program being stuck. This is the main reason I was thinking of implementing some kind of time out. I could also just fix my specific case in my program, but I thought a more general solution (which would probably look very similar anyways) to the main repo could be useful for everyone. I will try to dig into it later this week. If I find a (elegant) solution I will try to submit a pull request. I was thinking of using a timer signal from python's signal module, but it seems it only supports unix..

Either way, thanks for looking into it aswell, I really appreciate it!

[Jnelen]

After thinking about it some more, an alternative solution could involve having a parameter which allows you to set a limit for the maximum number of isomorphisms? If it is above the threshold, an error could be returned. The base case would be no limit, unless you think it'd be good to have a general ceiling? Either way, I still think a timeout feature would be a better solution, but the approach I mentioned here would be platform independent (and thus also work for windows), and should be relatively straight-forward to implement. What do you think @RMeli ? If you think it could be a good solution, I will give a try at implementing this later this week!

[RMeli]

an alternative solution could involve having a parameter which allows you to set a limit for the maximum number of isomorphisms?

Hello. This is what I was hinting at with the following comment:

I'm also wondering if there is some graph metric that can single out these pathological cases

However, I'm not aware if it's possible know a priori the number of isomorphisms (or to stop the networkx/graph-tool algorithm at a given number). Do you know of ways to compute (even approximately) the number of isomorphisms in a reasonable time?

I still think a timeout feature would be a better solution, but the approach I mentioned here would be platform independent

Yes, the timeout option is very interesting, but it would have to possibly be implemented for all platforms. One possible implementation would be using multiprocessing:

from multiprocessing import Pool, TimeoutError

from rdkit import Chem
from rdkit.Chem import AllChem

from spyrmsd.molecule import Molecule
from spyrmsd.rmsd import rmsdwrapper

# Benzene does not time out, RMSD is computed and printed
#rdmol1 = Chem.MolFromSmiles("c1ccccc1")
#rdmol2 = Chem.MolFromSmiles("c1ccccc1")

# Muparfostat times out
rdmol1 = Chem.MolFromSmiles("OP(O)(=O)OCC1OC(OC2C(OS(O)(=O)=O)C(COS(O)(=O)=O)OC(OC3C(OS(O)(=O)=O)OC(COS(O)(=O)=O)C(OS(O)(=O)=O)C3OS(O)(=O)=O)C2OS(O)(=O)=O)C(OS(O)(=O)=
O)C(OS(O)(=O)=O)C1OS(O)(=O)=O")
rdmol2 = Chem.MolFromSmiles("OP(O)(=O)OCC1OC(OC2C(OS(O)(=O)=O)C(COS(O)(=O)=O)OC(OC3C(OS(O)(=O)=O)OC(COS(O)(=O)=O)C(OS(O)(=O)=O)C3OS(O)(=O)=O)C2OS(O)(=O)=O)C(OS(O)(=O)=
O)C(OS(O)(=O)=O)C1OS(O)(=O)=O")

rdmol1 = AllChem.AddHs(rdmol1)
rdmol2 = AllChem.AddHs(rdmol2)
AllChem.EmbedMolecule(rdmol1)
AllChem.EmbedMolecule(rdmol2)

def compute_rmsd(rdmol1, rdmol2):
    mol1 = Molecule.from_rdkit(rdmol1)
    mol2 = Molecule.from_rdkit(rdmol2)

    return rmsdwrapper(mol1, mol2)
    

with Pool(processes=1) as pool:
    rmsd = pool.apply_async(compute_rmsd, (rdmol1, rdmol2))

    try:
        r = rmsd.get(timeout = 5)
        print(r)
    except TimeoutError:
        print("RMSD took too long to compute.")
        # Note: the task does not terminate
        # Note: the task terminates only when the pool context is exited

With benzene:

$ time python rmsd_mp.py
[0.5291976314926458]
python test.py  0.43s user 1.59s system 572% cpu 0.353 total

With muparfostat:

$ time python rmsd_mp.py
RMSD took too long to compute.
python test.py  4.31s user 3.27s system 128% cpu 5.913 total

However, with the overheads and additional complications that are involved here, I wonder if this should be implemented by users of spyrmsd instead of spyrmsd itself.

I encountered this problem before myself so it might be reasonably common, but it remains a bit of an edge case. What do you think?

[Jnelen]

Yes you're right. I (naively) assumed that determining the number of isomorphisms was a relatively quick process, thinking the real bottleneck would be calculating the RMSD for all the combinations. However, after some small-scale tests, it seems that determining the number of isomorphisms takes about half of the total compute time. Using that number to stop further calculations doesn't really solve the initial problem.
I think the multiprocessing solution is good in specific use-cases (like running on Windows or when you are executing things in parallel anyway), but as you mention it does create overhead which is probably not worth it in the vast majority of cases.
A good middle ground would have been to use signals, but again this is limited to Unix-based systems only. Since these really are edge-cases it's probably better to let users deal with it when necessary.

For my case, I'm probably going with the multiprocessing solution, as I was planning to add multicore support to my program anyway.
Thanks again for your responses and insights, they were really helpful and appreciated!

[Jnelen]

A short summary of solutions for future reference:

If you want to introduce some time-out functionality in your own code, there are 2 main options using native python modules: multiprocessing and signal.

• Multiprocessing
  Here you spawn a separate process for each calculation, where you can set a time-out.
  This works with Windows, and of course is a great option if you want to run things in parallel anyway, but it does introduce addtional overhead due to the spawning and managing of the processes. @RMeli provided a nice example earlier in this thread: Timeout functionality #105 (reply in thread)
• Signal
  Less overhead, but only available on Unix-based platforms. When the time runs out, a TimeoutError is raised, which you can handle accordingly to suit your needs. Code example (adapted from here):

from openbabel import pybel
from spyrmsd import io
from spyrmsd.molecule import Molecule
from spyrmsd.rmsd import rmsdwrapper

import signal
from contextlib import contextmanager

path1 = "DB01666_1.mol"
path2 = "DB01666_2.mol"

## Define TimeoutError
class TimeoutError(Exception):
    pass
  
@contextmanager
def timeout(time):

    def raise_timeout(signum, frame):
        raise TimeoutError("Ran out of time!")
    
    # Register a function to raise a TimeoutError on the signal.
    signal.signal(signal.SIGALRM, raise_timeout)
    # Schedule the signal to be sent after ``time``.
    signal.alarm(time)

    try:
        yield
    finally:
        # Unregister the signal so it won't be triggered
        # if the timeout is not reached.
        signal.signal(signal.SIGALRM, signal.SIG_IGN)
   
mol1 = io.loadmol(path1)
mol2 = io.loadmol(path2)

with timeout(5):
	rmsd = round(rmsdwrapper(mol1, mol2)[0],4)

print(rmsd)

[RMeli]

Thank you for the fantastic summary of our discussion @Jnelen, and for providing a MWE with signals.

I think the best solution for spyrmsd would be to have a high-level function where you can pass two list of molecules which computes the symmetry-corrected RMSD in parallel and allows to set a timeout (actually killing the working thread is needed in this case though, unlike the MWE I provided).

Unfortunately I don't have enough bandwidth to work on this.

For my case, I'm probably going with the multiprocessing solution, as I was planning to add multicore support to my program anyway.

If you end up implementing this solution, I'd be happy to have it integrated in spyrmsd as an high-level function (so that other people will not have to duplicate the same work). Some care is needed when using graph-tools, which uses OpenMP parallelisation behind the scenes.

[Jnelen]

Hi,

Thanks for pointing out the fact that graph-tools uses OpenMP, as I wasn't aware.
I noticed you can manually set the number of threads by using graph_tool.openmp_set_num_threads()
However at a first glance, i'm not really seeing a performance difference when setting it specifically to 1. Could it be that the functions spyrmsd uses aren't accelerated by multiple threads? I tried to run some of the "problematic" drugbank examples, and these were the results:

Using 1 cores for openMP
real    0m35.608s
user    0m34.465s
sys     0m3.927s

vs

Using 32 cores for openMP
real    0m35.110s
user    0m33.883s
sys     0m4.040s

Also, what would be the ideal input and output of the multiprocessing function? I was thinking a list of tuples/lists with length 2, which contain molecules/filepaths to the molecules? Other arguments would be number of cores(/processes) and an optional timeout? And for the output a list of RMSD values which maintains the order of the input lists?

I will play around with it and see what I can get working (well)!

[RMeli]

Ah. That's interesting. It's been a while since I looked into the graph isomorphism algorithm. It might be hard to parallelize properly. Incidentally, I was looking at rapids and noticed that they don't implement graph isomorphism either.

Also, what would be the ideal input and output of the multiprocessing function? I was thinking a list of tuples/lists with length 2, which contain molecules/filepaths to the molecules? Other arguments would be number of cores(/processes) and an optional timeout?

Yes, a list of tuples (or also two lists) with molecules would be ideal. And definitely with number of processes and (optional) timeout. Something like

def rmsdparallel(
    molrefs: List[molecule.Molecule],
    mols: List[molecule.Molecule],
    num_workers: int: 1,              # Maybe a better name for the argument?
    symmetry: bool = True,
    center: bool = False,
    minimize: bool = False,
    strip: bool = True,
    cache: bool = True,
    timeout: Optional[float] = None,
) -> Any:

that creates a multiprocessing pool and then dispatches to rmsdwrapper.

[RMeli]

I'm re-opening the discussion because GitHub only shows open discussions by default, and I don't want this to be hidden.

[Jnelen]

To follow-up on this, I have done quite a bit of experimenting. I ran into a lot of problems where I had difficulties finding a way to actually terminate the process (as you also indicated in your example earlier). However, I finally found something that seems to work decently well. If you think this is reasonable, I will polish it up and make a pull request. Basically I would add 2 functions to rmsd.py: one that just uses a timeout function on rmsdwrapper, and another one that is a multicore version of that. Ofcourse if you think it's more appropriate we could make the "regular" timeout function a private function as well. Anyway, here is the basic implementation:

def _compute_rmsd(mol, mol2, queue):
    rmsd = rmsdwrapper(mol1, mol2)[0]
    ## Use a Queue to store the result
    queue.put(rmsd)

def rmsd_with_timeout(sample, timeout=None):
        ## Inspired by https://superfastpython.com/task-with-timeout-child-process/
	mol1 = sample[0]
	mol2 = sample[1]
	
	queue = Queue()
	process = Process(target=_compute_rmsd, args=(mol1, mol2, queue))

	process.start()
	process.join(timeout=timeout)

	if process.is_alive():
		## terminate the process
		process.terminate()
		return None
	else:
		return queue.get()
		
def rmsd_timeout_multicore(samples, timeout=5, num_workers=1):    
	with ProcessPoolExecutor(max_workers=num_workers) as executor:
		rsmd_partial = partial(rmsd_with_timeout, timeout=timeout)
		results = executor.map(rsmd_partial,samples)
	
	return list(results)

I feel like this isn't too much code, and it only uses built-in python functions. From my small testing, the performance is still quite good (+ you can also use multiple cores now aswell). In the final implementation I would also include the other variables that rmsdwrapper can use, but this was to show the basic structure.

P.S.: naming suggestions are also welcome, I just put some indicative placeholders for now

[RMeli]

Thanks you very much for following up on this @Jnelen!

However, I finally found something that seems to work decently well. If you think this is reasonable, I will polish it up and make a pull request.

This sounds great, thank you for looking into this. I'll have to look a bit more closely into it, but from a first look it looks good.

Basically I would add 2 functions to rmsd.py: one that just uses a timeout function on rmsdwrapper, and another one that is a multicore version of that.

Do you think a single function would be enough? The multicore version could be used with n_workers=1 if one needs only the timeout and no parallelism.

TBH, on the long run, it might be worth having a single function that does everything, but this will require careful benchmarks and probably best to keep rmsdwrapper as it is for the time being.

I feel like this isn't too much code, and it only uses built-in python functions.

Yes, this looks quite clean/streamlined and does not require additional dependencies. I'd definitely be happy to have it within spyrmsd!

I'm quite bad at coming up with names... Maybe prmsdwrapper for parallel rmsdwrapper?

[Jnelen]

Great, thanks for the initial feedback! I opened a PR for this at #113. Let's start improving it from there!