Source Code for Module rdkit.Chem.FeatMaps.FeatMaps

  1  # $Id: FeatMaps.py 997 2009-02-25 06:12:43Z glandrum $ 
  2  # 
  3  # Copyright (C) 2006 Greg Landrum 
  4  # 
  5  #   @@ All Rights Reserved  @@ 
  6  # 
  7  from rdkit import Geometry 
  8  from rdkit import Chem 
  9  from rdkit.Chem import ChemicalFeatures 
 10  from rdkit.Chem.FeatMaps.FeatMapPoint import FeatMapPoint 
 11  import math 
 12   
 13 -class FeatMapScoreMode(object): 
 14    All=0 
 15    """ score each feature in the probe against every matching 
 16        feature in the FeatMap. 
 17    """ 
 18     
 19    Closest=1 
 20    """ score each feature in the probe against the closest 
 21        matching feature in the FeatMap. 
 22    """ 
 23   
 24    Best=2 
 25    """ score each feature in the probe against the matching 
 26        feature in the FeatMap that leads to the highest score 
 27    """ 
 28   
 29 -class FeatDirScoreMode(object): 
 30    Ignore=0 
 31    """ ignore feature directions 
 32    """ 
 33   
 34    DotFullRange=1 
 35    """ Use the dot product and allow negative contributions when 
 36        directions are anti-parallel. 
 37        e.g. score = dot(f1Dir,f2Dir) 
 38    """ 
 39   
 40    DotPosRange=2 
 41    """ Use the dot product and scale contributions to lie between 
 42        zero and one. 
 43        e.g. score = ( dot(f1Dir,f2Dir) + 1 ) / 2 
 44    """ 
 45   
 46 -class FeatMapParams(object): 
 47    """ one of these should be instantiated for each 
 48    feature type in the feature map 
 49    """ 
 50    radius=2.5 
 51    " cutoff radius " 
 52   
 53    width=1.0 
 54    " width parameter (e.g. the gaussian sigma) " 
 55     
 56 -  class FeatProfile(object): 
 57      " scoring profile of the feature " 
 58      Gaussian=0 
 59      Triangle=1 
 60      Box=2 
 61   
 62    featProfile=FeatProfile.Gaussian 
 63     
 64 -class FeatMap(object): 
 65    dirScoreMode = FeatDirScoreMode.Ignore 
 66    scoreMode=FeatMapScoreMode.All 
 67    params = {} 
 68 -  def __init__(self,params=None,feats=None, 
 69                 weights=None): 
 70      if params: 
 71        self.params = params 
 72      self._initializeFeats(feats,weights) 
 73       
 74 -  def _initializeFeats(self,feats,weights): 
 75      self._feats = [] 
 76      if feats: 
 77        if len(feats)!=len(weights): 
 78          raise ValueError,'feats and weights lists must be the same length' 
 79        for feat,weight in zip(feats,weights): 
 80          self.AddFeature(feat,weight) 
 81   
 82 -  def AddFeature(self,feat,weight=None): 
 83      if self.params and not self.params.has_key(feat.GetFamily()): 
 84        raise ValueError,'feature family %s not found in params'%feat.GetFamily() 
 85   
 86      newFeat = FeatMapPoint() 
 87      newFeat.initFromFeat(feat) 
 88      newFeat.weight = weight 
 89   
 90      self.AddFeatPoint(newFeat) 
 91       
 92 -  def AddFeatPoint(self,featPt): 
 93      if not isinstance(featPt,FeatMapPoint): 
 94        raise ValueError,'addFeatPoint() must be called with a FeatMapPoint instance' 
 95      if self.params and not self.params.has_key(featPt.GetFamily()): 
 96        raise ValueError,'feature family %s not found in params'%featPt.GetFamily() 
 97      self._feats.append(featPt) 
 98   
 99       
100 -  def GetFeatures(self): 
101      return self._feats 
102   
103 -  def GetNumFeatures(self): 
104      return len(self._feats) 
105   
106 -  def GetFeature(self, i): 
107      return self._feats[i] 
108   
109 -  def DropFeature(self,i): 
110      del self._feats[i] 
111   
112 -  def _loopOverMatchingFeats(self,oFeat): 
113      for sIdx,sFeat in enumerate(self._feats): 
114        if sFeat.GetFamily()==oFeat.GetFamily(): 
115          yield sIdx,sFeat 
116   
117 -  def GetFeatFeatScore(self,feat1,feat2,typeMatch=True): 
118      """ feat1 is one of our feats 
119          feat2 is any Feature 
120   
121   
122      """ 
123      if typeMatch and feat1.GetFamily()!=feat2.GetFamily(): 
124        return 0.0 
125      d2 = feat1.GetDist2(feat2) 
126      params = self.params[feat1.GetFamily()] 
127      if d2>params.radius*params.radius: 
128        return 0.0 
129   
130      if params.featProfile==FeatMapParams.FeatProfile.Gaussian: 
131        score = math.exp(-d2/params.width) 
132      elif params.featProfile==FeatMapParams.FeatProfile.Triangle: 
133        d = math.sqrt(d2) 
134        if d<params.width: 
135          score = 1.-d/params.width 
136        else: 
137          score=0.0 
138      elif params.featProfile==FeatMapParams.FeatProfile.Box: 
139        score = 1.0 
140      weight = feat1.weight 
141      score *= weight 
142   
143      if self.dirScoreMode != FeatDirScoreMode.Ignore: 
144        dirScore = feat1.GetDirMatch(feat2) 
145        if self.dirScoreMode==FeatDirScoreMode.DotPosRange: 
146          dirScore = (dirScore + 1.0)/2.0 
147        elif self.dirScoreMode!=FeatDirScoreMode.DotFullRange: 
148          raise NotImplementedError,'bad feature dir score mode' 
149        score *= dirScore 
150         
151      return score 
152       
153 -  def ScoreFeats(self,featsToScore,mapScoreVect=[],featsScoreVect=[], 
154                   featsToFeatMapIdx=[]): 
155      nFeats = len(self._feats) 
156      if mapScoreVect and len(mapScoreVect)!=nFeats: 
157        raise ValueError,'if provided, len(mapScoreVect) should equal numFeats' 
158      nToScore = len(featsToScore) 
159      if featsScoreVect and len(featsScoreVect)!=nToScore: 
160        raise ValueError,'if provided, len(featsScoreVect) should equal len(featsToScore)' 
161      if featsToFeatMapIdx and len(featsToFeatMapIdx)!=nToScore: 
162        raise ValueError,'if provided, len(featsToFeatMapIdx) should equal len(featsToScore)' 
163   
164      if mapScoreVect: 
165        for i in range(nFeats): mapScoreVect[i]=0.0 
166      else: 
167        mapScoreVect = [0.0]*nFeats 
168         
169      if self.scoreMode==FeatMapScoreMode.Closest: 
170        defScore=1000.0 
171      else: 
172        defScore=0.0 
173      if featsScoreVect: 
174        for i in range(nToScore): featsScoreVect[i]=defScore 
175      else: 
176        featsScoreVect = [defScore]*nToScore 
177         
178      if not featsToFeatMapIdx: 
179        featsToFeatMapIdx = [None]*nToScore 
180        
181      for i in range(nToScore): 
182        if self.scoreMode != FeatMapScoreMode.All: 
183          featsToFeatMapIdx[i]=[-1] 
184        else: 
185          featsToFeatMapIdx[i]=[] 
186         
187      for oIdx,oFeat in enumerate(featsToScore): 
188        for sIdx,sFeat in self._loopOverMatchingFeats(oFeat): 
189          if self.scoreMode == FeatMapScoreMode.Closest: 
190            d = sFeat.GetDist2(oFeat) 
191            if d<featsScoreVect[oIdx]: 
192              featsScoreVect[oIdx]=d 
193              featsToFeatMapIdx[oIdx][0]=sIdx 
194          else: 
195            lScore = self.GetFeatFeatScore(sFeat,oFeat,typeMatch=False) 
196            if self.scoreMode == FeatMapScoreMode.Best: 
197              if lScore>featsScoreVect[oIdx]: 
198                featsScoreVect[oIdx]=lScore 
199                featsToFeatMapIdx[oIdx][0]=sIdx 
200            elif self.scoreMode == FeatMapScoreMode.All: 
201              featsScoreVect[oIdx] += lScore 
202              mapScoreVect[sIdx] += lScore 
203              featsToFeatMapIdx[oIdx].append(sIdx) 
204            else: 
205              raise ValueError,'bad score mode' 
206   
207      totScore = 0.0 
208      if self.scoreMode == FeatMapScoreMode.Closest: 
209        for oIdx,oFeat in enumerate(featsToScore): 
210          sIdx = featsToFeatMapIdx[oIdx][0] 
211          if sIdx>-1: 
212            lScore = self.GetFeatFeatScore(sFeat,oFeat,typeMatch=False) 
213            featsScoreVect[oIdx] = lScore 
214            mapScoreVect[sIdx] = lScore 
215            totScore += lScore 
216          else: 
217            featsScoreVect[oIdx] = 0 
218             
219      else: 
220        totScore = sum(featsScoreVect) 
221        if self.scoreMode == FeatMapScoreMode.Best: 
222          for oIdx,lScore in enumerate(featsScoreVect): 
223            sIdx = featsToFeatMapIdx[oIdx][0] 
224            if sIdx>-1: 
225              mapScoreVect[sIdx]=lScore 
226   
227      # replace placeholders: 
228      if self.scoreMode != FeatMapScoreMode.All: 
229        for elem in featsToFeatMapIdx: 
230          if elem==[-1]: 
231            elem.pop() 
232      return totScore 
233       
234   
235 -  def __str__(self): 
236      res = '' 
237      for i,feat in enumerate(self._feats): 
238        weight = feat.weight 
239        pos = feat.GetPos() 
240        res += '% 3d % 12s % 6.4f % 6.4f % 6.4f % 6.4f\n'%(i+1, 
241                                                           feat.GetFamily(), 
242                                                           pos.x,pos.y,pos.z,weight) 
243      return res 
244   
245   
246  #------------------------------------ 
247  # 
248  #  doctest boilerplate 
249  # 
250 -def _test(): 
251    import doctest,sys 
252    return doctest.testmod(sys.modules["__main__"]) 
253   
254  if __name__ == '__main__': 
255    import sys 
256    failed,tried = _test() 
257    sys.exit(failed) 
258