Package rdkit :: Package Chem :: Package Subshape :: Module BuilderUtils
[hide private]
[frames] | no frames]

Source Code for Module rdkit.Chem.Subshape.BuilderUtils

  1  # 
  2  # Copyright (C) 2007-2017 by Greg Landrum 
  3  #  All rights reserved 
  4  # 
  5  from __future__ import print_function 
  6   
  7  import math 
  8   
  9  import numpy 
 10   
 11  from rdkit import Geometry 
 12  from rdkit.Chem.Subshape import SubshapeObjects 
 13   
 14   
15 -def ComputeGridIndices(shapeGrid, winRad):
16 if getattr(shapeGrid, '_indicesInSphere', None): 17 return shapeGrid._indicesInSphere 18 gridSpacing = shapeGrid.GetSpacing() 19 dX = shapeGrid.GetNumX() 20 dY = shapeGrid.GetNumY() 21 radInGrid = int(winRad / gridSpacing) 22 indicesInSphere = [] 23 for i in range(-radInGrid, radInGrid + 1): 24 for j in range(-radInGrid, radInGrid + 1): 25 for k in range(-radInGrid, radInGrid + 1): 26 d = int(math.sqrt(i * i + j * j + k * k)) 27 if d <= radInGrid: 28 idx = (i * dY + j) * dX + k 29 indicesInSphere.append(idx) 30 shapeGrid._indicesInSphere = indicesInSphere 31 return indicesInSphere
32 33
34 -def ComputeShapeGridCentroid(pt, shapeGrid, winRad):
35 count = 0 36 centroid = Geometry.Point3D(0, 0, 0) 37 idxI = shapeGrid.GetGridPointIndex(pt) 38 shapeGridVect = shapeGrid.GetOccupancyVect() 39 40 indicesInSphere = ComputeGridIndices(shapeGrid, winRad) 41 42 nGridPts = len(shapeGridVect) 43 for idxJ in indicesInSphere: 44 idx = idxI + idxJ 45 if idx >= 0 and idx < nGridPts: 46 wt = shapeGridVect[idx] 47 tPt = shapeGrid.GetGridPointLoc(idx) 48 centroid += tPt * wt 49 count += wt 50 if not count: 51 raise ValueError('found no weight in sphere') 52 centroid /= count 53 # print 'csgc:','(%2f,%2f,%2f)'%tuple(pt),'(%2f,%2f,%2f)'%tuple(centroid),count 54 return count, centroid
55 56
57 -def FindTerminalPtsFromShape(shape, winRad, fraction, maxGridVal=3):
58 pts = Geometry.FindGridTerminalPoints(shape.grid, winRad, fraction) 59 termPts = [SubshapeObjects.SkeletonPoint(location=x) for x in pts] 60 return termPts
61 62
63 -def FindTerminalPtsFromConformer(conf, winRad, nbrCount):
64 mol = conf.GetOwningMol() 65 nAts = conf.GetNumAtoms() 66 nbrLists = [[] for _ in range(nAts)] 67 for i in range(nAts): 68 if (mol.GetAtomWithIdx(i).GetAtomicNum() <= 1): 69 continue 70 pi = conf.GetAtomPosition(i) 71 nbrLists[i].append((i, pi)) 72 for j in range(i + 1, nAts): 73 if (mol.GetAtomWithIdx(j).GetAtomicNum() <= 1): 74 continue 75 pj = conf.GetAtomPosition(j) 76 dist = pi.Distance(conf.GetAtomPosition(j)) 77 if dist < winRad: 78 nbrLists[i].append((j, pj)) 79 nbrLists[j].append((i, pi)) 80 termPts = [] 81 # for i in range(nAts): 82 # if not len(nbrLists[i]): continue 83 # if len(nbrLists[i])>10: 84 # print i+1,len(nbrLists[i]) 85 # else: 86 # print i+1,len(nbrLists[i]),[x[0]+1 for x in nbrLists[i]] 87 88 while 1: 89 for i in range(nAts): 90 if not nbrLists[i]: 91 continue 92 pos = Geometry.Point3D(0, 0, 0) 93 totWt = 0.0 94 if len(nbrLists[i]) < nbrCount: 95 nbrList = nbrLists[i] 96 for j in range(0, len(nbrList)): 97 nbrJ, posJ = nbrList[j] 98 weight = 1. * len(nbrLists[i]) / len(nbrLists[nbrJ]) 99 pos += posJ * weight 100 totWt += weight 101 pos /= totWt 102 termPts.append(SubshapeObjects.SkeletonPoint(location=pos)) 103 if not len(termPts): 104 nbrCount += 1 105 else: 106 break 107 return termPts
108 109
110 -def FindGridPointBetweenPoints(pt1, pt2, shapeGrid, winRad):
111 center = pt1 + pt2 112 center /= 2.0 113 d = 1e8 114 while d > shapeGrid.GetSpacing(): 115 count, centroid = Geometry.ComputeGridCentroid(shapeGrid, center, winRad) 116 d = center.Distance(centroid) 117 center = centroid 118 return center
119 120
121 -def ClusterTerminalPts(pts, winRad, scale):
122 res = [] 123 tagged = [(y, x) for x, y in enumerate(pts)] 124 while tagged: 125 head, headIdx = tagged.pop(0) 126 currSet = [head] 127 i = 0 128 while i < len(tagged): 129 nbr, nbrIdx = tagged[i] 130 if head.location.Distance(nbr.location) < scale * winRad: 131 currSet.append(nbr) 132 del tagged[i] 133 else: 134 i += 1 135 pt = Geometry.Point3D(0, 0, 0) 136 for o in currSet: 137 pt += o.location 138 pt /= len(currSet) 139 res.append(SubshapeObjects.SkeletonPoint(location=pt)) 140 return res
141 142
143 -def GetMoreTerminalPoints(shape, pts, winRad, maxGridVal, targetNumber=5):
144 """ adds a set of new terminal points using a max-min algorithm 145 """ 146 shapeGrid = shape.grid 147 shapeVect = shapeGrid.GetOccupancyVect() 148 nGridPts = len(shapeVect) 149 # loop, taking the grid point with the maximum minimum distance, until 150 # we have enough points 151 while len(pts) < targetNumber: 152 maxMin = -1 153 for i in range(nGridPts): 154 if shapeVect[i] < maxGridVal: 155 continue 156 minVal = 1e8 157 posI = shapeGrid.GetGridPointLoc(i) 158 for currPt in pts: 159 dst = posI.Distance(currPt.location) 160 if dst < minVal: 161 minVal = dst 162 if minVal > maxMin: 163 maxMin = minVal 164 bestPt = posI 165 count, centroid = Geometry.ComputeGridCentroid(shapeGrid, bestPt, winRad) 166 pts.append(SubshapeObjects.SkeletonPoint(location=centroid))
167 168
169 -def FindFarthestGridPoint(shape, loc, winRad, maxGridVal):
170 """ find the grid point with max occupancy that is furthest from a 171 given location 172 """ 173 shapeGrid = shape.grid 174 shapeVect = shapeGrid.GetOccupancyVect() 175 nGridPts = len(shapeVect) 176 dMax = -1 177 for i in range(nGridPts): 178 if shapeVect[i] < maxGridVal: 179 continue 180 posI = shapeGrid.GetGridPointLoc(i) 181 dst = posI.Distance(loc) 182 if dst > dMax: 183 dMax = dst 184 res = posI 185 186 count, centroid = Geometry.ComputeGridCentroid(shapeGrid, res, winRad) 187 res = centroid 188 return res
189 190
191 -def ExpandTerminalPts(shape, pts, winRad, maxGridVal=3.0, targetNumPts=5):
192 """ find additional terminal points until a target number is reached 193 """ 194 if len(pts) == 1: 195 # if there's only one point, find the grid point with max value that is 196 # *farthest* from this one and use it: 197 pt2 = FindFarthestGridPoint(shape, pts[0].location, winRad, maxGridVal) 198 pts.append(SubshapeObjects.SkeletonPoint(location=pt2)) 199 if len(pts) == 2: 200 # add a point roughly in the middle: 201 shapeGrid = shape.grid 202 pt1 = pts[0].location 203 pt2 = pts[1].location 204 center = FindGridPointBetweenPoints(pt1, pt2, shapeGrid, winRad) 205 pts.append(SubshapeObjects.SkeletonPoint(location=center)) 206 if len(pts) < targetNumPts: 207 GetMoreTerminalPoints(shape, pts, winRad, maxGridVal, targetNumPts)
208 209
210 -def AppendSkeletonPoints(shapeGrid, termPts, winRad, stepDist, maxGridVal=3, maxDistC=15.0, 211 distTol=1.5, symFactor=1.5, verbose=False):
212 nTermPts = len(termPts) 213 skelPts = [] 214 shapeVect = shapeGrid.GetOccupancyVect() 215 nGridPts = len(shapeVect) 216 # find all possible skeleton points 217 if verbose: 218 print('generate all possible') 219 for i in range(nGridPts): 220 if shapeVect[i] < maxGridVal: 221 continue 222 posI = shapeGrid.GetGridPointLoc(i) 223 ok = True 224 for pt in termPts: 225 dst = posI.Distance(pt.location) 226 if dst < stepDist: 227 ok = False 228 break 229 if ok: 230 skelPts.append(SubshapeObjects.SkeletonPoint(location=posI)) 231 # now start removing them 232 if verbose: 233 print('Compute centroids:', len(skelPts)) 234 gridBoxVolume = shapeGrid.GetSpacing()**3 235 maxVol = 4.0 * math.pi / 3.0 * winRad**3 * maxGridVal / gridBoxVolume 236 i = 0 237 while i < len(skelPts): 238 pt = skelPts[i] 239 count, centroid = Geometry.ComputeGridCentroid(shapeGrid, pt.location, winRad) 240 # count,centroid=ComputeShapeGridCentroid(pt.location,shapeGrid,winRad) 241 centroidPtDist = centroid.Distance(pt.location) 242 if centroidPtDist > maxDistC: 243 del skelPts[i] 244 else: 245 pt.fracVol = float(count) / maxVol 246 pt.location.x = centroid.x 247 pt.location.y = centroid.y 248 pt.location.z = centroid.z 249 i += 1 250 251 if verbose: 252 print('remove points:', len(skelPts)) 253 res = termPts + skelPts 254 i = 0 255 while i < len(res): 256 p = -1 257 mFrac = 0.0 258 ptI = res[i] 259 260 startJ = max(i + 1, nTermPts) 261 for j in range(startJ, len(res)): 262 ptJ = res[j] 263 distC = ptI.location.Distance(ptJ.location) 264 if distC < symFactor * stepDist: 265 if ptJ.fracVol > mFrac: 266 p = j 267 mFrac = ptJ.fracVol 268 # print i,len(res),p,mFrac 269 if p > -1: 270 ptP = res.pop(p) 271 j = startJ 272 while j < len(res): 273 ptJ = res[j] 274 distC = ptI.location.Distance(ptJ.location) 275 if distC < symFactor * stepDist: 276 del res[j] 277 else: 278 j += 1 279 res.append(ptP) 280 # print '% 3d'%i,'% 5.2f % 5.2f % 5.2f'%tuple(list(ptI.location)),' - ','% 5.2f % 5.2f % 5.2f'%tuple(list(ptJ.location)) 281 i += 1 282 return res
283 284
285 -def CalculateDirectionsAtPoint(pt, shapeGrid, winRad):
286 shapeGridVect = shapeGrid.GetOccupancyVect() 287 nGridPts = len(shapeGridVect) 288 tmp = winRad / shapeGrid.GetSpacing() 289 radInGrid = int(tmp) 290 radInGrid2 = int(tmp * tmp) 291 covMat = numpy.zeros((3, 3), numpy.float64) 292 293 dX = shapeGrid.GetNumX() 294 dY = shapeGrid.GetNumY() 295 # dZ = shapeGrid.GetNumZ() 296 idx = shapeGrid.GetGridPointIndex(pt.location) 297 idxZ = idx // (dX * dY) 298 rem = idx % (dX * dY) 299 idxY = rem // dX 300 idxX = rem % dX 301 totWt = 0.0 302 for i in range(-radInGrid, radInGrid): 303 xi = idxX + i 304 for j in range(-radInGrid, radInGrid): 305 xj = idxY + j 306 for k in range(-radInGrid, radInGrid): 307 xk = idxZ + k 308 d2 = i * i + j * j + k * k 309 if d2 > radInGrid2 and int(math.sqrt(d2)) > radInGrid: 310 continue 311 gridIdx = (xk * dY + xj) * dX + xi 312 if gridIdx >= 0 and gridIdx < nGridPts: 313 wtHere = shapeGridVect[gridIdx] 314 totWt += wtHere 315 ptInSphere = shapeGrid.GetGridPointLoc(gridIdx) 316 ptInSphere -= pt.location 317 covMat[0][0] += wtHere * (ptInSphere.x * ptInSphere.x) 318 covMat[0][1] += wtHere * (ptInSphere.x * ptInSphere.y) 319 covMat[0][2] += wtHere * (ptInSphere.x * ptInSphere.z) 320 covMat[1][1] += wtHere * (ptInSphere.y * ptInSphere.y) 321 covMat[1][2] += wtHere * (ptInSphere.y * ptInSphere.z) 322 covMat[2][2] += wtHere * (ptInSphere.z * ptInSphere.z) 323 covMat /= totWt 324 covMat[1][0] = covMat[0][1] 325 covMat[2][0] = covMat[0][2] 326 covMat[2][1] = covMat[1][2] 327 328 eVals, eVects = numpy.linalg.eigh(covMat) 329 sv = list(zip(eVals, numpy.transpose(eVects))) 330 sv.sort(reverse=True) 331 eVals, eVects = list(zip(*sv)) 332 pt.shapeMoments = tuple(eVals) 333 pt.shapeDirs = tuple([Geometry.Point3D(p[0], p[1], p[2]) for p in eVects])
334 335 # print '-------------' 336 # print pt.location.x,pt.location.y,pt.location.z 337 # for v in covMat: 338 # print ' ',v 339 # print '---' 340 # print eVals 341 # for v in eVects: 342 # print ' ',v 343 # print '-------------' 344 345
346 -def AssignMolFeatsToPoints(pts, mol, featFactory, winRad):
347 feats = featFactory.GetFeaturesForMol(mol) 348 for i, pt in enumerate(pts): 349 for feat in feats: 350 if feat.GetPos().Distance(pt.location) < winRad: 351 typ = feat.GetFamily() 352 if typ not in pt.molFeatures: 353 pt.molFeatures.append(typ) 354 print(i, pt.molFeatures)
355