Reference implementation for computing myopic Maximum Common Edge Subgraph (MCES) distances, described in the publication Coverage bias in small molecule machine learning (Citation).
Note
New in versions >= 2.1.0: much faster computations with the C++ implementation of the lower bounds (default bound ~266× faster than Python, see benchmarks). Install via pip install myopic-mces[cpp].
Input and Output file are in csv format per default. Every line in the input-file is one comparison:
input-file: index,SMILES1,SMILES2
output-file: index,myopic MCES distance,computation time in seconds,computation mode
The output file contains additional information if the MCES mapping for the structure and/or numbered smiles are computed:
output-file: index,myopic MCES distance,computation time in seconds,computation mode,mapping
output-file: index,myopic MCES distance,computation time in seconds,computation mode,mapping,numsmiles1,numsmiles2
output-file: index,myopic MCES distance,computation time in seconds,computation mode,numsmiles1,numsmiles2
If the MCES structure is computed, it is saved as a mapping. If the numsmiles are saved as well, they contain the indexed smiles in which each atom is numbered. The structure mapping has the format {(AtomID_SMILES1, AtomID_SMILES1): (AtomID_SMILES2, AtomID_SMILES2)}, where the first AtomID pair is a bond in the molecule depicted by SMILES1 and the second pair is a bond in the molecule depicted by
SMILES2.
Download via pip and execute:
pip install myopic-mces
myopic_mces input-file [output-file]
The `output-file` argument may be omitted when using `--hdf5_mode` (see below).The lower bounds are also available as a much faster C++ implementation, which is used automatically when installed (see Recommended settings):
pip install myopic-mces[cpp]Alternatively, to install directly from this repository:
pip install -e .For usage in Python:
from myopic_mces import MCES
MCES('CC(=O)OC1=CC=CC=C1C(=O)O', 'CN1C=NC2=C1C(=O)N(C(=O)N2C)C')See the PuLP documentation on how to configure ILP solvers.
As PuLP does not ship with solvers anymore, a solver must be configured before running myopic-mces.
The default solver for myopic-mces is 'COIN_CMD', which can be installed via:
pip install pulp[cbc]
General options
--threshold int Threshold for the comparison.
Exact distance is only calculated if the distance is lower than the threshold.
If set to -1 the exact distance is always calculated.
Default: 10
--solver string Solver used for solving the ILP. Examples:'CPLEX_CMD', 'GUROBI_CMD', 'GLPK_CMD'.
Default: COIN_CMD
--num_jobs int Number of jobs; instances to run in parallel.
By default this is set to the number of (logical) CPU cores.
Default: -1
--hdf5_mode If set, input will be read from `input-file` in HDF5 format; output will be appended to this file.
In this mode the `output-file` positional argument is not needed and may be omitted.
See "Prepare HDF5 input" below
--hide_rdkit_warnings If set, attempts to supress RDKit warnings.
Options for the ILP solver
--solver_onethreaded If set, limits ILP solver to one thread, likely resulting in faster
performance with parallel computations (not available for all solvers).
--solver_no_msg If set, prevents solver from logging (not available for all solvers)
--solver_time_limit_seconds Set a time limit for the ILP solver. Instances that hit the time limit get the
computation mode 5 (solution found, not proven exact) or 6 (no solution found,
distance -1). Only supported for the solvers COIN_CMD and CPLEX_PY.
With --solver_time_limit_seconds (or timeLimit in solver_options), the computation mode of each instance tells whether its distance can be trusted:
| mode | meaning |
|---|---|
| 1 | exact distance, proven optimal by the solver |
| 2 | distance is above the threshold, the threshold is returned |
| 5 | time limit hit, a solution was found but not proven exact (upper bound of the distance) |
| 6 | time limit hit, no solution found, distance -1 |
This distinction is only reliable for the solvers COIN_CMD and CPLEX_PY, so a time limit is rejected for all other solvers. For example, CPLEX_CMD reports a time limit without solution as infeasible, which would wrongly be returned as mode 2 (above threshold). To use a time limit, switch to one of the supported solvers.
Similarly, solver gap tolerances that are too loose are rejected (gapAbs must be below 0.5, gapRel must not be above 1e-4),
as solutions within these tolerances would otherwise be reported as exact (mode 1).
Experimental options for myopic MCES distance computation
--no_ilp_threshold If set, do not add threshold as constraint to ILP,
resulting in longer runtimes and potential violations of the triangle equation.
--choose_bound_dynamically If set, a potentially weaker but faster lower bound will be computed and used
when this bound is already greater than the threshold. By default (without
this option), always the strongest lower bound will be computed and used.
--use_bound_zero If set, compute and use an additional weaker formula-based lower bound.
Use in conjunction with `choose_bound_dynamically`.
--catch_computation_errors Instead of aborting the computation, instances that failed to compute receive distance "-1".
--jobs_batch_size Batch size for parallelization.
Default: 32
--jobs_dispatch Pre-dispatch of jobs for parallelization.
Default: 10*n_jobs
--use_matrix_lookup Use with the path to a HDF5 file with precomputed MCES distances. Computation for these
instances will be skipped, using the provided values. HDF5 has to contain distances (key `mces`)
and SMILES (`mces_smiles_order`), like the HDF5 files produced by this script.
NOTE: When used in combination with `prepare_input`, only use with `--no_shuffle`.
--lookup_threshold Use with `--use_matrix_lookup`: Precomputed values equal or greater than the threshold
will be ignored; these instances will be recomputed.
--num_smiles If set, calculates and saves SMILES of numbered atoms to the output for each input smiles.
--structure If set, saves a mapping of the pairwise MCES structure to the output. Only calculated for molecules
for which the MCES distance has been proven exact by the ILP.
To speed up computations and save space, use the variant of the package implemented in C++ (install via pip install myopic_mces[cpp]; computes the default lower bound ~266× faster than the Python implementation, see benchmarks), the CPLEX solver (requires pip install myopic_mces[cplex] for newer CPLEX versions), HDF5-mode (see below) and enable --choose_bound_dynamically:
PATH=$CPLEX_HOME/bin/x86-64_linux/:$PATH python -m myopic_mces.myopic_mces --threshold 10 --solver CPLEX_CMD --solver_onethreaded --solver_no_msg --hdf5_mode input-file.hdf5The PATH-variable might have to be adapted, depending on the CPLEX version, to contain the directory of the CPLEX executable (see the PuLP documentation).
Python packages required are:
rdkit(==2022.09.5)
networkx(==3.0)
pulp[cbc](==2.7.0)
scipy(==1.10.1)
joblib(==1.2.0)
Version numbers in braces correspond to an exemplary tested configuration (under Python version 3.11.0). The program can be run on any standard operating system, tested on Windows 10 64 bit and Arch-Linux@linux-6.2.7 64 bit.
The recommended method of installation is via pip.
pip install myopic_mcesOptional extras enable additional functionality:
pip install myopic_mces[hdf5] # for HDF5 mode
pip install myopic_mces[cplex] # for newer versions of the CPLEX solver
pip install myopic_mces[cpp] # for the fast C++ implementation of the lower boundsImportant
The C++ bounds are built against a specific RDKit version (currently rdkit==2025.3.4), and installing the cpp extra pins RDKit to that version. Other RDKit versions are not supported at the moment. If you cannot use this RDKit version, use the pure Python implementation (install without the cpp extra).
When working from a clone of this repository (e.g., for development), the recommended method is uv. Download this repository, navigate to the download location and execute the following commands:
uv sync
# optionally with extras, e.g.:
uv sync --extra hdf5 --extra cpp
# to activate the created environment:
source .venv/bin/activateAlternatively, conda/mamba can be used, although pip or uv are recommended.
The provided environment includes the hdf5 and cpp extras (RDKit is installed via pip, as required by the C++ bounds):
conda env create -f conda_env.yml # or: mamba env create -f conda_env.yml
conda activate myopic_mces
pip install cplex # optional, for the CPLEX solverA typical installation time should not exceed 5 minutes, mostly depending on the internet connection speed to download all required packages.
The example provided in example/example_data.csv can be run with:
pip install myopic-mces
myopic_mces example/example_data.csv example/example_data_out.csvTypical runtime is about 10s on Windows 10 with all default options, running on 4 cores with 8GB RAM. Exemplary output is provided in example/example_data_out.csv.
For big datasets it is recommended to divide the input into batches, which can be done with this script (HDF5 mode recommended, see below!):
python -m myopic_mces.prepare_input input-file.csv output-folder/ --batch_size 50_000_000input-file.csv has to be formatted as shown above. This creates output_folder with the subdirectory data puts all batches (named batch$i.csv) inside.
Requires pip install myopic_mces[hdf5].
To conserve space, input for myopic MCES computation can now be provided as a HDF5-file containing the following "Datasets":
smiles: list of all unique SMILEScomputation_indices: matrix with the shape (n, 3), each row representing one instance to compute. The first column contains the index of the pair, the second and third the indices of the two SMILES, respectively
Example:
smiles = [smiles_a, smiles_b, smiles_c]
computation_indices = [[0, 0, 1], # computes MCES for a vs. b
[1, 0, 2], # computes MCES for a vs. c
[2, 1, 2], # computes MCES for b vs. c
]
Instead of preparing the HDF5-file manually, with the additional option to create batches, prepare_input.py can be used in HDF5-mode:
python -m myopic_mces.prepare_input --batch_size 50_000_000 --hdf5_mode input-smiles.txt output-folder/
# or for two datasets, dataset_a versus dataset_b:
python -m myopic_mces.prepare_input --batch_size 50_000_000 --hdf5_mode --hdf5_extra_input_file dataset-b-smiles.txt dataset-a-smiles.txt output-folder/Created batches (batch$i.hdf5) are written directly to the output-folder.
With --use_db_lookup, pairs already known in a database (mcesdb, not public yet) are reused and land in batch0_precomputed.hdf5; batches to compute start with batch1.hdf5. Database connection parameters come from the PGHOST/PGDATABASE/... environment variables.
python -m myopic_mces.prepare_input --batch_size 50_000_000 --hdf5_mode --use_db_lookup --threshold 10 input-smiles.txt output-folder/Requires pip install myopic_mces[hdf5].
Batches results can be combined into a single HDF5 file with this script:
python -m myopic_mces.combine_hdf5_batches --out output-file batches/batch*.hdf5To get a square matrix of MCES distances for convenience within Python from this file, simply do:
import h5py
from scipy.spatial.distance import squareform
with h5py.File(hdf5outputfile, 'r') as f:
mces_square = squareform(f['mces'])
smiles = f['mces_smiles_order'] # column/row labelsIf you want to filter datasets for similar structures in a database, you can use this script. This can speed up computations considerably, as for each query structure computations are stopped, when a "match" is found.
python filter_dataset.py --input input-file.json --out output-file --threshold 10Input is provided in the following format (JSON file):
{query_smiles1: [db_smiles1, db_smiles2, ...],
query_smiles2: [db_smiles1, db_smiles2, ...]}
The database can be different for each query SMILES, allow pre-filtering depending on the query.
F. Kretschmer, J. Seipp, M. Ludwig, G. W. Klau, and S. Böcker. Coverage bias in small molecule machine learning. Nat Commun 16(1):554, 2025.