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

Source Code for Module rdkit.Chem.Subshape.BuilderUtils

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