Package rdkit :: Package Chem :: Package FeatMaps :: Module FeatMaps
[hide private]
[frames] | no frames]

Source Code for Module rdkit.Chem.FeatMaps.FeatMaps

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