Chem.Subshape.BuilderUtils

1 # $Id: BuilderUtils.py 285 2007-07-15 15:00:36Z glandrum $ 2 # 3 # Copyright (C) 2007 by Greg Landrum 4 # All rights reserved 5 # 6 import Geometry 7 from Chem.Subshape import SubshapeObjects 8 import math 9 import Numeric 10 import LinearAlgebra 11 12 13 #-----------------------------------------------------------------------------

14 -def ComputeGridIndices(shapeGrid,winRad):

15 if getattr(shapeGrid,'_indicesInSphere',None): 16 return shapeGrid._indicesInSphere 17 gridSpacing = shapeGrid.GetSpacing() 18 dX = shapeGrid.GetNumX() 19 dY = shapeGrid.GetNumY() 20 dZ = shapeGrid.GetNumZ() 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 58 #-----------------------------------------------------------------------------

59 -def FindTerminalPtsFromShape(shape,winRad,fraction,maxGridVal=3):

60 pts = Geometry.FindGridTerminalPoints(shape.grid,winRad,fraction) 61 termPts = [SubshapeObjects.SkeletonPoint(location=x) for x in pts] 62 return termPts

63 64 #-----------------------------------------------------------------------------

65 -def FindTerminalPtsFromConformer(conf,winRad,nbrCount):

66 mol = conf.GetOwningMol() 67 nAts = conf.GetNumAtoms() 68 nbrLists=[[] for x in range(nAts)] 69 for i in range(nAts): 70 if(mol.GetAtomWithIdx(i).GetAtomicNum()<=1): continue 71 pi = conf.GetAtomPosition(i) 72 nbrLists[i].append((i,pi)) 73 for j in range(i+1,nAts): 74 if(mol.GetAtomWithIdx(j).GetAtomicNum()<=1): 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]: continue 91 pos = Geometry.Point3D(0,0,0) 92 totWt=0.0 93 if len(nbrLists[i])<nbrCount: 94 nbrList = nbrLists[i] 95 for j in range(0,len(nbrList)): 96 nbrJ,posJ=nbrList[j] 97 weight = 1.*len(nbrLists[i])/len(nbrLists[nbrJ]) 98 pos += posJ*weight 99 totWt+=weight 100 pos /= totWt 101 termPts.append(SubshapeObjects.SkeletonPoint(location=pos)) 102 if not len(termPts): 103 nbrCount += 1 104 else: 105 break 106 return termPts

107 108 #-----------------------------------------------------------------------------

109 -def FindGridPointBetweenPoints(pt1,pt2,shapeGrid,winRad):

110 center = pt1+pt2 111 center /= 2.0 112 d=1e8 113 while d>shapeGrid.GetSpacing(): 114 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,center,winRad) 115 d = center.Distance(centroid) 116 center = centroid 117 return center

118 119 #-----------------------------------------------------------------------------

120 -def ClusterTerminalPts(pts,winRad,scale):

121 res = [] 122 tagged = [(y,x) for x,y in enumerate(pts)] 123 while tagged: 124 head,headIdx = tagged.pop(0) 125 currSet = [head] 126 i=0 127 while i<len(tagged): 128 nbr,nbrIdx=tagged[i] 129 if head.location.Distance(nbr.location)<scale*winRad: 130 currSet.append(nbr) 131 del tagged[i] 132 else: 133 i+=1 134 pt = Geometry.Point3D(0,0,0) 135 for o in currSet: 136 pt += o.location 137 pt /= len(currSet) 138 res.append(SubshapeObjects.SkeletonPoint(location=pt)) 139 return res

140

141 -def GetMoreTerminalPoints(shape,pts,winRad,maxGridVal,targetNumber=5):

142 """ adds a set of new terminal points using a max-min algorithm 143 """ 144 shapeGrid=shape.grid 145 shapeVect = shapeGrid.GetOccupancyVect() 146 nGridPts = len(shapeVect) 147 # loop, taking the grid point with the maximum minimum distance, until 148 # we have enough points 149 while len(pts)<targetNumber: 150 maxMin=-1 151 for i in range(nGridPts): 152 if shapeVect[i]<maxGridVal: 153 continue 154 minVal=1e8 155 posI = shapeGrid.GetGridPointLoc(i) 156 for currPt in pts: 157 dst = posI.Distance(currPt.location) 158 if dst<minVal: 159 minVal=dst 160 if minVal>maxMin: 161 maxMin=minVal 162 bestPt=posI 163 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,bestPt,winRad) 164 pts.append(SubshapeObjects.SkeletonPoint(location=centroid))

165 166

167 -def FindFarthestGridPoint(shape,loc,winRad,maxGridVal):

168 """ find the grid point with max occupancy that is furthest from a 169 given location 170 """ 171 shapeGrid=shape.grid 172 shapeVect = shapeGrid.GetOccupancyVect() 173 nGridPts = len(shapeVect) 174 dMax=-1; 175 for i in range(nGridPts): 176 if shapeVect[i]<maxGridVal: 177 continue 178 posI = shapeGrid.GetGridPointLoc(i) 179 dst = posI.Distance(loc) 180 if dst>dMax: 181 dMax=dst 182 res=posI 183 184 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,res,winRad) 185 res=centroid 186 return res

187

188 -def ExpandTerminalPts(shape,pts,winRad,maxGridVal=3.0,targetNumPts=5):

189 """ find additional terminal points until a target number is reached 190 """ 191 if len(pts)==1: 192 # if there's only one point, find the grid point with max value that is 193 # *farthest* from this one and use it: 194 pt2=FindFarthestGridPoint(shape,pts[0].location,winRad,maxGridVal) 195 pts.append(SubshapeObjects.SkeletonPoint(location=pt2)) 196 if len(pts)==2: 197 # add a point roughly in the middle: 198 shapeGrid=shape.grid 199 pt1 = pts[0].location 200 pt2 = pts[1].location 201 center = FindGridPointBetweenPoints(pt1,pt2,shapeGrid,winRad) 202 pts.append(SubshapeObjects.SkeletonPoint(location=center)) 203 if len(pts)<targetNumPts: 204 GetMoreTerminalPoints(shape,pts,winRad,maxGridVal,targetNumPts)

205 206 207 #-----------------------------------------------------------------------------

208 -def AppendSkeletonPoints(shapeGrid,termPts,winRad,stepDist,maxGridVal=3, 209 maxDistC=15.0,distTol=1.5,symFactor=1.5):

210 nTermPts = len(termPts) 211 skelPts=[] 212 shapeVect = shapeGrid.GetOccupancyVect() 213 nGridPts = len(shapeVect) 214 # find all possible skeleton points 215 print 'generate all possible' 216 for i in range(nGridPts): 217 if shapeVect[i]<maxGridVal: 218 continue 219 posI = shapeGrid.GetGridPointLoc(i) 220 ok=True 221 for pt in termPts: 222 dst = posI.Distance(pt.location) 223 if dst<stepDist: 224 ok=False 225 break 226 if ok: 227 skelPts.append(SubshapeObjects.SkeletonPoint(location=posI)) 228 # now start removing them 229 print 'Compute centroids:',len(skelPts) 230 gridBoxVolume=shapeGrid.GetSpacing()**3 231 maxVol = 4.0*math.pi/3.0 * winRad**3 * maxGridVal / gridBoxVolume 232 i=0 233 while i<len(skelPts): 234 pt = skelPts[i] 235 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,pt.location,winRad) 236 #count,centroid=ComputeShapeGridCentroid(pt.location,shapeGrid,winRad) 237 centroidPtDist=centroid.Distance(pt.location) 238 if centroidPtDist>maxDistC: 239 del skelPts[i] 240 else: 241 pt.fracVol = float(count)/maxVol 242 pt.location.x = centroid.x 243 pt.location.y = centroid.y 244 pt.location.z = centroid.z 245 i+=1 246 247 print 'remove points:',len(skelPts) 248 res = termPts+skelPts 249 i=0 250 while i<len(res): 251 p=-1 252 mFrac=0.0 253 ptI = res[i] 254 255 startJ = max(i+1,nTermPts) 256 for j in range(startJ,len(res)): 257 ptJ=res[j] 258 distC = ptI.location.Distance(ptJ.location) 259 if distC<symFactor*stepDist: 260 if ptJ.fracVol>mFrac: 261 p=j 262 mFrac=ptJ.fracVol 263 #print i,len(res),p,mFrac 264 if p>-1: 265 ptP = res.pop(p) 266 j = startJ 267 while j < len(res): 268 ptJ=res[j] 269 distC = ptI.location.Distance(ptJ.location) 270 if distC<symFactor*stepDist: 271 del res[j] 272 else: 273 j+=1 274 res.append(ptP) 275 #print '% 3d'%i,'% 5.2f % 5.2f % 5.2f'%tuple(list(ptI.location)),' - ','% 5.2f % 5.2f % 5.2f'%tuple(list(ptJ.location)) 276 i+=1 277 return res

278 279 #-----------------------------------------------------------------------------

280 -def CalculateDirectionsAtPoint(pt,shapeGrid,winRad):

281 shapeGridVect = shapeGrid.GetOccupancyVect() 282 nGridPts = len(shapeGridVect) 283 tmp = winRad/shapeGrid.GetSpacing() 284 radInGrid=int(tmp) 285 radInGrid2=int(tmp*tmp) 286 covMat = Numeric.zeros((3,3),Numeric.Float) 287 288 dX = shapeGrid.GetNumX() 289 dY = shapeGrid.GetNumY() 290 dZ = shapeGrid.GetNumZ() 291 idx = shapeGrid.GetGridPointIndex(pt.location) 292 idxZ = idx//(dX*dY) 293 rem = idx%(dX*dY) 294 idxY = rem//dX 295 idxX = rem%dX 296 totWt=0.0 297 for i in range(-radInGrid,radInGrid): 298 xi = idxX+i 299 for j in range(-radInGrid,radInGrid): 300 xj = idxY+j 301 for k in range(-radInGrid,radInGrid): 302 xk = idxZ+k 303 d2 = i*i+j*j+k*k 304 if d2>radInGrid2 and int(math.sqrt(d2))>radInGrid: 305 continue 306 gridIdx = (xk*dY+xj)*dX+xi 307 if gridIdx>=0 and gridIdx<nGridPts: 308 wtHere = shapeGridVect[gridIdx] 309 totWt += wtHere 310 ptInSphere = shapeGrid.GetGridPointLoc(gridIdx) 311 ptInSphere -= pt.location 312 covMat[0][0]+= wtHere*(ptInSphere.x*ptInSphere.x) 313 covMat[0][1]+= wtHere*(ptInSphere.x*ptInSphere.y) 314 covMat[0][2]+= wtHere*(ptInSphere.x*ptInSphere.z) 315 covMat[1][1]+= wtHere*(ptInSphere.y*ptInSphere.y) 316 covMat[1][2]+= wtHere*(ptInSphere.y*ptInSphere.z) 317 covMat[2][2]+= wtHere*(ptInSphere.z*ptInSphere.z) 318 covMat /= totWt 319 covMat[1][0] = covMat[0][1] 320 covMat[2][0] = covMat[0][2] 321 covMat[2][1] = covMat[1][2] 322 323 eVals,eVects = LinearAlgebra.eigenvectors(covMat) 324 sv = [(x,y) for x,y in zip(eVals,eVects)] 325 sv.sort(reverse=True) 326 eVals = [x for x,y in sv] 327 eVects = [y for x,y in 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

Source Code for Module Chem.Subshape.BuilderUtils