1
2
3
4
5
6
7
8
9
10
11 """ utility functionality for the 2D pharmacophores code
12
13 See Docs/Chem/Pharm2D.triangles.jpg for an illustration of the way
14 pharmacophores are broken into triangles and labelled.
15
16 See Docs/Chem/Pharm2D.signatures.jpg for an illustration of bit
17 numbering
18
19 """
20 from __future__ import print_function, division
21 import numpy
22
23
24
25
26
27 nPointDistDict = {
28 2: ((0,1),),
29 3: ((0,1),(0,2),
30 (1,2)),
31 4: ((0,1),(0,2),(0,3),
32 (1,2),(2,3)),
33 5: ((0,1),(0,2),(0,3),(0,4),
34 (1,2),(2,3),(3,4)),
35 6: ((0,1),(0,2),(0,3),(0,4),(0,5),
36 (1,2),(2,3),(3,4),(4,5)),
37 7: ((0,1),(0,2),(0,3),(0,4),(0,5),(0,6),
38 (1,2),(2,3),(3,4),(4,5),(5,6)),
39 8: ((0,1),(0,2),(0,3),(0,4),(0,5),(0,6),(0,7),
40 (1,2),(2,3),(3,4),(4,5),(5,6),(6,7)),
41 9: ((0,1),(0,2),(0,3),(0,4),(0,5),(0,6),(0,7),(0,8),
42 (1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8)),
43 10: ((0,1),(0,2),(0,3),(0,4),(0,5),(0,6),(0,7),(0,8),(0,9),
44 (1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9)),
45 }
46
47
48
49
50 nDistPointDict = {
51 1:2,
52 3:3,
53 5:4,
54 7:5,
55 9:6,
56 11:7,
57 13:8,
58 15:9,
59 17:10,
60 }
61
62 _trianglesInPharmacophore = {}
81
82
84 if x <= 1:
85 return 1
86
87 accum = 1
88 for i in range(x):
89 accum *= i+1
90 return accum
91
93 """ checks the triangle inequality for combinations
94 of distance bins.
95
96 the general triangle inequality is:
97 d1 + d2 >= d3
98 the conservative binned form of this is:
99 d1(upper) + d2(upper) >= d3(lower)
100
101 """
102 if d1[1]+d2[1]<d3[0]: return False
103 if d2[1]+d3[1]<d1[0]: return False
104 if d3[1]+d1[1]<d2[0]: return False
105
106 return True
107
109 """ checks the scaffold passed in to see if all
110 contributing triangles can satisfy the triangle inequality
111
112 the scaffold itself (encoded in combo) is a list of binned distances
113
114 """
115
116 nPts = nDistPointDict[len(combo)]
117 tris = GetTriangles(nPts)
118 for tri in tris:
119 ds = [bins[combo[x]] for x in tri]
120 if not BinsTriangleInequality(ds[0],ds[1],ds[2]):
121 return False
122 return True
123
124
125 _numCombDict = {}
127 """ returns the number of ways to fit nItems into nSlots
128
129 We assume that (x,y) and (y,x) are equivalent, and
130 (x,x) is allowed.
131
132 General formula is, for N items and S slots:
133 res = (N+S-1)! / ( (N-1)! * S! )
134
135 """
136 global _numCombDict
137 res = _numCombDict.get((nItems,nSlots),-1)
138 if res == -1:
139 res = _fact(nItems+nSlots-1) // (_fact(nItems-1)*_fact(nSlots))
140 _numCombDict[(nItems,nSlots)] = res
141 return res
142
143 _verbose = 0
144
145 _countCache={}
146 -def CountUpTo(nItems,nSlots,vs,idx=0,startAt=0):
147 """ Figures out where a given combination of indices would
148 occur in the combinatorial explosion generated by _GetIndexCombinations_
149
150 **Arguments**
151
152 - nItems: the number of items to distribute
153
154 - nSlots: the number of slots in which to distribute them
155
156 - vs: a sequence containing the values to find
157
158 - idx: used in the recursion
159
160 - startAt: used in the recursion
161
162 **Returns**
163
164 an integer
165
166 """
167 global _countCache
168 if _verbose:
169 print(' '*idx,'CountUpTo(%d)'%idx,vs[idx],startAt)
170 if idx==0 and (nItems,nSlots,tuple(vs)) in _countCache:
171 return _countCache[(nItems,nSlots,tuple(vs))]
172 elif idx >= nSlots:
173 accum = 0
174 elif idx == nSlots-1:
175 accum = vs[idx]-startAt
176 else:
177 accum = 0
178
179 for i in range(startAt,vs[idx]):
180 nLevsUnder = nSlots-idx-1
181 nValsOver = nItems-i
182 if _verbose:
183 print(' '*idx,' ',i,nValsOver,nLevsUnder,
184 NumCombinations(nValsOver,nLevsUnder))
185 accum += NumCombinations(nValsOver,nLevsUnder)
186 accum += CountUpTo(nItems,nSlots,vs,idx+1,vs[idx])
187 if _verbose: print(' '*idx,'>',accum)
188 if idx == 0:
189 _countCache[(nItems,nSlots,tuple(vs))] = accum
190 return accum
191
192 _indexCombinations={}
194 """ Generates all combinations of nItems in nSlots without including
195 duplicates
196
197 **Arguments**
198
199 - nItems: the number of items to distribute
200
201 - nSlots: the number of slots in which to distribute them
202
203 - slot: used in recursion
204
205 - lastItemVal: used in recursion
206
207 **Returns**
208
209 a list of lists
210
211 """
212 global _indexCombinations
213 if not slot and (nItems,nSlots) in _indexCombinations:
214 res = _indexCombinations[(nItems,nSlots)]
215 elif slot >= nSlots:
216 res = []
217 elif slot == nSlots-1:
218 res = [[x] for x in range(lastItemVal,nItems)]
219 else:
220 res = []
221 for x in range(lastItemVal,nItems):
222 tmp = GetIndexCombinations(nItems,nSlots,slot+1,x)
223 for entry in tmp:
224 res.append([x]+entry)
225 if not slot:
226 _indexCombinations[(nItems,nSlots)] = res
227 return res
228
230 """ Does the combinatorial explosion of the possible combinations
231 of the elements of _choices_.
232
233 **Arguments**
234
235 - choices: sequence of sequences with the elements to be enumerated
236
237 - noDups: (optional) if this is nonzero, results with duplicates,
238 e.g. (1,1,0), will not be generated
239
240 - which: used in recursion
241
242 **Returns**
243
244 a list of lists
245
246 >>> GetAllCombinations([(0,),(1,),(2,)])
247 [[0, 1, 2]]
248 >>> GetAllCombinations([(0,),(1,3),(2,)])
249 [[0, 1, 2], [0, 3, 2]]
250
251 >>> GetAllCombinations([(0,1),(1,3),(2,)])
252 [[0, 1, 2], [0, 3, 2], [1, 3, 2]]
253
254 """
255 if which >= len(choices):
256 res = []
257 elif which == len(choices)-1:
258 res = [[x] for x in choices[which]]
259 else:
260 res = []
261 tmp = GetAllCombinations(choices,noDups=noDups,
262 which=which+1)
263 for thing in choices[which]:
264 for other in tmp:
265 if not noDups or thing not in other:
266 res.append([thing]+other)
267 return res
268
270 """ Does the combinatorial explosion of the possible combinations
271 of the elements of _choices_.
272
273 """
274 assert len(choices)==len(classes)
275 if which >= len(choices):
276 res = []
277 elif which == len(choices)-1:
278 res = [[(classes[which],x)] for x in choices[which]]
279 else:
280 res = []
281 tmp = GetUniqueCombinations(choices,classes,
282 which=which+1)
283 for thing in choices[which]:
284 for other in tmp:
285 idxThere=0
286 for x in other:
287 if x[1]==thing:idxThere+=1
288 if not idxThere:
289 newL = [(classes[which],thing)]+other
290 newL.sort()
291 if newL not in res:
292 res.append(newL)
293 return res
294
296 """ uniquifies the combinations in the argument
297
298 **Arguments**:
299
300 - combos: a sequence of sequences
301
302 **Returns**
303
304 - a list of tuples containing the unique combos
305
306 """
307 print('>>> u:',combos)
308 resD = {}
309 for combo in combos:
310 k = combo[:]
311 k.sort()
312 resD[tuple(k)] = tuple(combo)
313 print(' >>> u:',resD.values())
314 return resD.values()
315
317 """ gets all realizable scaffolds (passing the triangle inequality) with the
318 given number of points and returns them as a list of tuples
319
320 """
321 if nPts < 2:
322 res = 0
323 elif nPts == 2:
324 res = [(x,) for x in range(len(bins))]
325 else:
326 nDists = len(nPointDistDict[nPts])
327 combos = GetAllCombinations([range(len(bins))]*nDists,noDups=0)
328 res = []
329 for combo in combos:
330 if not useTriangleInequality or ScaffoldPasses(combo,bins):
331 res.append(tuple(combo))
332 return res
333
335 """
336 put the distances for a triangle into canonical order
337
338 It's easy if the features are all different:
339 >>> OrderTriangle([0,2,4],[1,2,3])
340 ([0, 2, 4], [1, 2, 3])
341
342 It's trickiest if they are all the same:
343 >>> OrderTriangle([0,0,0],[1,2,3])
344 ([0, 0, 0], [3, 2, 1])
345 >>> OrderTriangle([0,0,0],[2,1,3])
346 ([0, 0, 0], [3, 2, 1])
347 >>> OrderTriangle([0,0,0],[1,3,2])
348 ([0, 0, 0], [3, 2, 1])
349 >>> OrderTriangle([0,0,0],[3,1,2])
350 ([0, 0, 0], [3, 2, 1])
351 >>> OrderTriangle([0,0,0],[3,2,1])
352 ([0, 0, 0], [3, 2, 1])
353
354 >>> OrderTriangle([0,0,1],[3,2,1])
355 ([0, 0, 1], [3, 2, 1])
356 >>> OrderTriangle([0,0,1],[1,3,2])
357 ([0, 0, 1], [1, 3, 2])
358 >>> OrderTriangle([0,0,1],[1,2,3])
359 ([0, 0, 1], [1, 3, 2])
360 >>> OrderTriangle([0,0,1],[1,3,2])
361 ([0, 0, 1], [1, 3, 2])
362
363 """
364 if len(featIndices)!=3: raise ValueError('bad indices')
365 if len(dists)!=3: raise ValueError('bad dists')
366
367 fs = set(featIndices)
368 if len(fs)==3:
369 return featIndices,dists
370
371 dSums=[0]*3
372 dSums[0] = dists[0]+dists[1]
373 dSums[1] = dists[0]+dists[2]
374 dSums[2] = dists[1]+dists[2]
375 mD = max(dSums)
376 if len(fs)==1:
377 if dSums[0]==mD:
378 if dists[0]>dists[1]:
379 ireorder=(0,1,2)
380 dreorder=(0,1,2)
381 else:
382 ireorder=(0,2,1)
383 dreorder=(1,0,2)
384 elif dSums[1]==mD:
385 if dists[0]>dists[2]:
386 ireorder=(1,0,2)
387 dreorder=(0,2,1)
388 else:
389 ireorder=(1,2,0)
390 dreorder=(2,0,1)
391 else:
392 if dists[1]>dists[2]:
393 ireorder = (2,0,1)
394 dreorder=(1,2,0)
395 else:
396 ireorder = (2,1,0)
397 dreorder=(2,1,0)
398 else:
399
400 if featIndices[0]==featIndices[1]:
401 if dists[1]>dists[2]:
402 ireorder=(0,1,2)
403 dreorder=(0,1,2)
404 else:
405 ireorder=(1,0,2)
406 dreorder=(0,2,1)
407 elif featIndices[0]==featIndices[2]:
408 if dists[0]>dists[2]:
409 ireorder=(0,1,2)
410 dreorder=(0,1,2)
411 else:
412 ireorder=(2,1,0)
413 dreorder=(2,1,0)
414 else:
415 if dists[0]>dists[1]:
416 ireorder=(0,1,2)
417 dreorder=(0,1,2)
418 else:
419 ireorder=(0,2,1)
420 dreorder=(1,0,2)
421 dists = [dists[x] for x in dreorder]
422 featIndices = [featIndices[x] for x in ireorder]
423 return featIndices,dists
424
425
426
427
428
430 import doctest,sys
431 return doctest.testmod(sys.modules["__main__"])
432
433
434 if __name__ == '__main__':
435 import sys
436 failed,tried = _test()
437 sys.exit(failed)
438