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