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