1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
//
// Copyright (C) 2020 Greg Landrum and T5 Informatics GmbH
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include "MolEnumerator.h"
#include <RDGeneral/Exceptions.h>
namespace RDKit {
namespace MolEnumerator {
namespace {
//! recursively builds the variations
void getVariations(size_t level, std::vector<size_t> base,
std::vector<std::vector<size_t>> &variations,
const std::vector<size_t> &variationCounts,
size_t maxToEnumerate, bool doRandom) {
PRECONDITION(level < variationCounts.size(), "bad recursion");
for (size_t i = 0; i < variationCounts[level]; ++i) {
base[level] = i;
if (level + 1 >= variationCounts.size()) {
// at the bottom of the recursion
variations.push_back(base);
} else {
getVariations(level + 1, base, variations, variationCounts,
maxToEnumerate, doRandom);
}
if (variations.size() >= maxToEnumerate) {
return;
}
}
}
void enumerateVariations(std::vector<std::vector<size_t>> &variations,
const std::vector<size_t> &variationCounts,
const MolEnumeratorParams ¶ms) {
if (params.doRandom) {
UNDER_CONSTRUCTION("random enumeration not yet supported");
}
variations.clear();
std::vector<size_t> base(variationCounts.size(), 0);
getVariations(0, base, variations, variationCounts, params.maxToEnumerate,
params.doRandom);
}
} // namespace
namespace detail {
const std::string idxPropName = "_enumeratorOrigIdx";
void preserveOrigIndices(ROMol &mol) {
mol.setProp(idxPropName, 1);
for (auto atom : mol.atoms()) {
atom->setProp(idxPropName, atom->getIdx());
}
}
void removeOrigIndices(ROMol &mol) {
for (auto atom : mol.atoms()) {
atom->clearProp(idxPropName);
}
mol.clearProp(idxPropName);
}
} // namespace detail
namespace {
void clearReactionProps(ROMol &mol) {
bool includeRings = false;
mol.clearComputedProps(includeRings);
for (auto atom : mol.atoms()) {
atom->clearProp(common_properties::reactantAtomIdx);
atom->clearProp(common_properties::reactantIdx);
atom->clearProp("was_dummy");
atom->clearProp(common_properties::reactionMapNum);
}
}
} // namespace
MolBundle enumerate(const ROMol &mol,
const std::vector<MolEnumeratorParams> ¶mLists) {
std::unique_ptr<MolBundle> accum{new MolBundle()};
boost::shared_ptr<ROMol> molCpy{new ROMol(mol)};
detail::preserveOrigIndices(*molCpy);
accum->addMol(molCpy);
bool variationsFound = false;
for (const auto ¶ms : paramLists) {
if (!params.dp_operation) {
throw ValueErrorException("MolEnumeratorParams has no operation set");
}
std::unique_ptr<MolBundle> thisRound{new MolBundle()};
for (const auto &tmol : accum->getMols()) {
// copy the op since we will modify it:
auto op = params.dp_operation->copy();
op->initFromMol(*tmol);
auto variationCounts = op->getVariationCounts();
if (variationCounts.empty()) {
RDKit::ROMOL_SPTR mcp(new ROMol(*tmol));
clearReactionProps(*mcp);
thisRound->addMol(mcp);
continue;
}
std::vector<std::vector<size_t>> variations;
enumerateVariations(variations, variationCounts, params);
for (const auto &variation : variations) {
variationsFound = true;
auto newMol = (*op)(variation);
newMol->updatePropertyCache(false);
clearReactionProps(*newMol);
thisRound->addMol(ROMOL_SPTR(newMol.release()));
}
}
accum.swap(thisRound);
}
if (!variationsFound) {
return MolBundle();
}
for (auto rmol : accum->getMols()) {
detail::removeOrigIndices(*rmol);
}
return *accum;
}
MolBundle enumerate(const ROMol &mol, size_t maxPerOperation) {
std::vector<MolEnumeratorParams> paramsList;
// position variation first
MolEnumerator::MolEnumeratorParams posVariationParams;
auto pvOp = new MolEnumerator::PositionVariationOp();
posVariationParams.dp_operation =
std::shared_ptr<MolEnumerator::MolEnumeratorOp>(pvOp);
if (maxPerOperation > 0) {
posVariationParams.maxToEnumerate = maxPerOperation;
}
paramsList.push_back(posVariationParams);
// now repeat units
MolEnumerator::MolEnumeratorParams repeatUnitParams;
auto sruOp = new MolEnumerator::RepeatUnitOp();
repeatUnitParams.dp_operation =
std::shared_ptr<MolEnumerator::MolEnumeratorOp>(sruOp);
if (maxPerOperation > 0) {
sruOp->d_maxNumRounds = maxPerOperation;
repeatUnitParams.maxToEnumerate = maxPerOperation;
}
paramsList.push_back(repeatUnitParams);
// linknodes last because we can only enumerate mols with a single
// fragment
MolEnumerator::MolEnumeratorParams linkParams;
auto lpOp = new MolEnumerator::LinkNodeOp();
linkParams.dp_operation =
std::shared_ptr<MolEnumerator::MolEnumeratorOp>(lpOp);
if (maxPerOperation > 0) {
linkParams.maxToEnumerate = maxPerOperation;
}
paramsList.push_back(linkParams);
return enumerate(mol, paramsList);
}
} // namespace MolEnumerator
} // namespace RDKit
|
