# $Id$
#
#  Copyright (C) 2005-2006 Rational Discovery LLC
#
#   @@ All Rights Reserved @@
#  This file is part of the RDKit.
#  The contents are covered by the terms of the BSD license
#  which is included in the file license.txt, found at the root
#  of the RDKit source tree.
#
_version = "$Rev$"
_splashMessage="""
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  FeatFinderCLI version %s

  Copyright (C) 2005 Rational Discovery LLC

  This software is copyrighted.  The software may not be copied,
  reproduced, translated or reduced to any electronic medium or
  machine-readable form without the prior written consent of
  Rational Discovery LLC.
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
"""%_version
from rdkit import Chem
from rdkit.Chem import ChemicalFeatures
from rdkit import RDLogger
logger = RDLogger.logger()
import sys,os
import re
splitExpr = re.compile(r'[ \t,]')

def GetAtomFeatInfo(factory,mol):
  res = [None]*mol.GetNumAtoms()
  feats = factory.GetFeaturesForMol(mol)
  for feat in feats:
    ids = feat.GetAtomIds()
    for id in ids:
      if res[id] is None:
        res[id] = []
      res[id].append("%s-%s"%(feat.GetFamily(),feat.GetType()))
  return res

if __name__ == '__main__':
  def Usage():
    message="""
    Usage: FeatFinderCLI [-r] <fdefFilename> <smilesFilename>
    
    NOTE:
      - the smiles file should have SMILES in the first column

    """
    print >>sys.stderr,message


  import getopt
  args,extras = getopt.getopt(sys.argv[1:],'r')
  reverseIt=False
  for arg,val in args:
    if arg=='-r':
      reverseIt=True
      
  if len(extras)<2:
    Usage()
    sys.exit(-1)
  print >>sys.stderr,_splashMessage
  fdefFilename = extras[0]
  if not os.path.exists(fdefFilename):
    logger.error("Fdef file %s does not exist."%fdefFilename)
    sys.exit(-1)
  try:
    factory = ChemicalFeatures.BuildFeatureFactory(fdefFilename)
  except:
    logger.error("Could not parse Fdef file %s."%fdefFilename,exc_info=True)
    sys.exit(-1)
    
  smilesFilename = extras[1]
  if not os.path.exists(smilesFilename):
    logger.error("Smiles file %s does not exist."%smilesFilename)
    sys.exit(-1)

  try:
    inF = file(smilesFilename,'r')
  except:
    logger.error("Could not open smiles file %s."%smilesFilename,exc_info=True)
    sys.exit(-1)

  lineNo=0
  for line in inF.readlines():
    lineNo+=1
    line = line.strip()
    smi = splitExpr.split(line)[0].strip()
    try:
      mol = Chem.MolFromSmiles(smi)
    except:
      mol = None

    if mol:
      print 'Mol-%d\t%s'%(lineNo,smi)

      if not reverseIt:
        featInfo = GetAtomFeatInfo(factory,mol)
        for i,v in enumerate(featInfo):
          print '\t% 2s(%d)'%(mol.GetAtomWithIdx(i).GetSymbol(),i+1),
          if v:
            print '\t',', '.join(v)
          else:
            print
      else:
        feats = factory.GetFeaturesForMol(mol)
        for feat in feats:
          print '\t%s-%s: '%(feat.GetFamily(),feat.GetType()),
          print ', '.join([str(x) for x in feat.GetAtomIds()])
    else:
      logger.warning("Could not process smiles '%s' on line %d."%(smi,lineNo))