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