// $Id$
//
//  Copyright (C) 2001-2010 Greg Landrum and 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.
//
#include <GraphMol/RDKitBase.h>
#include <GraphMol/Canon.h>
#include <RDGeneral/Exceptions.h>
#include <RDGeneral/hash/hash.hpp>
#include <algorithm>

namespace RDKit {
namespace Canon {
  struct _possibleCompare : public std::binary_function<PossibleType,PossibleType,bool> {
    bool operator()(const PossibleType &arg1,const PossibleType &arg2) const {
      return (arg1.get<0>() < arg2.get<0>());
    }
  };

  void switchBondDir(Bond *bond){
    PRECONDITION(bond,"bad bond");
    PRECONDITION(bond->getBondType()==Bond::SINGLE||bond->getIsAromatic(),
                 "bad bond type");
    switch(bond->getBondDir()){
    case Bond::ENDUPRIGHT:
      bond->setBondDir(Bond::ENDDOWNRIGHT);
      break;
    case Bond::ENDDOWNRIGHT:
      bond->setBondDir(Bond::ENDUPRIGHT);
      break;
    default:
      break;
    }
  }
  
  // FIX: this may only be of interest from the SmilesWriter, should we
  // move it there?
  //
  //
  void canonicalizeDoubleBond(Bond *dblBond,
                              INT_VECT &bondVisitOrders,
                              INT_VECT &atomVisitOrders,
                              INT_VECT &bondDirCounts,
                              INT_VECT &atomDirCounts){
    PRECONDITION(dblBond,"bad bond");
    PRECONDITION(dblBond->getBondType() == Bond::DOUBLE,"bad bond order");
    PRECONDITION(dblBond->getStereo() > Bond::STEREOANY,"bad bond stereo");
    PRECONDITION(dblBond->getStereoAtoms().size() >= 2,"bad bond stereo atoms");
    PRECONDITION(atomVisitOrders[dblBond->getBeginAtomIdx()]>0 ||
                 atomVisitOrders[dblBond->getEndAtomIdx()]>0,
                 "neither end atom traversed");
    
    // atom1 is the lower numbered atom of the double bond (the one traversed
    // first)
    Atom *atom1,*atom2;
    if(atomVisitOrders[dblBond->getBeginAtomIdx()] <
       atomVisitOrders[dblBond->getEndAtomIdx()] ){
      atom1 = dblBond->getBeginAtom();
      atom2 = dblBond->getEndAtom();
    } else {
      atom1 = dblBond->getEndAtom();
      atom2 = dblBond->getBeginAtom();
    }

    // we only worry about double bonds that begin and end at atoms
    // of degree 2 or 3:
    if( (atom1->getDegree() != 2 && atom1->getDegree() != 3) ||
        (atom2->getDegree() != 2 && atom2->getDegree() != 3) ) {
      return;
    }
    
    Bond *firstFromAtom1=NULL,*secondFromAtom1=NULL;
    Bond *firstFromAtom2=NULL,*secondFromAtom2=NULL;

    int firstVisitOrder=100000;

    ROMol &mol=dblBond->getOwningMol();
    
    ROMol::OBOND_ITER_PAIR atomBonds;
    // -------------------------------------------------------
    // find the lowest visit order bonds from each end and determine
    // if anything is already constraining our choice of directions:
    bool dir1Set=false,dir2Set=false;
    atomBonds = mol.getAtomBonds(atom1);
    while( atomBonds.first != atomBonds.second ){
      if(mol[*atomBonds.first].get() != dblBond){
        int bondIdx = mol[*atomBonds.first]->getIdx();      
        if( bondDirCounts[bondIdx] > 0 ){
          dir1Set = true;
        }
        if(!firstFromAtom1 || bondVisitOrders[bondIdx] < firstVisitOrder ){
          if(firstFromAtom1) secondFromAtom1 = firstFromAtom1;
          firstFromAtom1 = mol[*atomBonds.first].get();
          firstVisitOrder = bondVisitOrders[bondIdx];
        } else {
          secondFromAtom1 = mol[*atomBonds.first].get();
        }
      }
      atomBonds.first++;
    }
    atomBonds = mol.getAtomBonds(atom2);
    firstVisitOrder = 10000;
    while( atomBonds.first != atomBonds.second ){
      if(mol[*atomBonds.first].get() != dblBond){
        int bondIdx =mol[*atomBonds.first]->getIdx();      
        if( bondDirCounts[bondIdx] > 0 ){
          dir2Set = true;
        }
        if(!firstFromAtom2 || bondVisitOrders[bondIdx] < firstVisitOrder){
          if(firstFromAtom2) secondFromAtom2 = firstFromAtom2;
          firstFromAtom2 = mol[*atomBonds.first].get();
          firstVisitOrder = bondVisitOrders[bondIdx];
        } else {
          secondFromAtom2 = mol[*atomBonds.first].get();
        }
      }
      atomBonds.first++;
    }

    // make sure we found everything we need to find:
    CHECK_INVARIANT(firstFromAtom1,"could not find atom1");
    CHECK_INVARIANT(firstFromAtom2,"could not find atom2");
    CHECK_INVARIANT(atom1->getDegree()==2 || secondFromAtom1,"inconsistency at atom1");
    CHECK_INVARIANT(atom2->getDegree()==2 || secondFromAtom2,"inconsistency at atom2");

    bool setFromBond1=true;
    Bond::BondDir atom1Dir=Bond::NONE;
    Bond::BondDir atom2Dir=Bond::NONE;
    Bond *atom1ControllingBond=firstFromAtom1;
    Bond *atom2ControllingBond=firstFromAtom2;
    if(!dir1Set && !dir2Set){
      // ----------------------------------
      // nothing has touched our bonds so far, so set the
      // directions to "arbitrary" values:

      // the bond we came in on becomes ENDUPRIGHT:
      atom1Dir = Bond::ENDUPRIGHT;
      firstFromAtom1->setBondDir(atom1Dir);
      bondDirCounts[firstFromAtom1->getIdx()] += 1;
      atomDirCounts[atom1->getIdx()]+=1;
    } else if(!dir2Set){
      // at least one of the bonds on atom1 has its directionality set already:
      if(bondDirCounts[firstFromAtom1->getIdx()]>0){
        // The first bond's direction has been set at some earlier point:
        atom1Dir = firstFromAtom1->getBondDir();
        bondDirCounts[firstFromAtom1->getIdx()] += 1;
        atomDirCounts[atom1->getIdx()]+=1;
        if(secondFromAtom1){
          // both bonds have their directionalities set, make sure
          // they are compatible:
          if( firstFromAtom1->getBondDir()==secondFromAtom1->getBondDir() &&
              bondDirCounts[firstFromAtom2->getIdx()] ){
            CHECK_INVARIANT(((firstFromAtom1->getBeginAtomIdx()==atom1->getIdx()) ^
                             (secondFromAtom1->getBeginAtomIdx()==atom1->getIdx()))
                            ,"inconsistent state");
          }
        }
      } else {
        // the second bond must be present and setting the direction:
        CHECK_INVARIANT(secondFromAtom1,"inconsistent state");
        CHECK_INVARIANT(bondDirCounts[secondFromAtom1->getIdx()]>0,"inconsistent state");
        // It must be the second bond setting the direction.
        // This happens when the bond dir is set in a branch:
        //        v- this double bond
        //   CC(/C=P/N)=N/O
        //      ^- the second bond sets the direction
        // or when the first bond is a ring closure from an
        // earlier traversed atom:
        //             v- this double bond
        //   NC1=NOC/C1=N\O
        //     ^- this closure ends up being the first bond,
        //        and it does not set the direction.
        //
        // This addresses parts of Issue 185 and sf.net Issue 1842174
        // 
        atom1Dir = secondFromAtom1->getBondDir();

        firstFromAtom1->setBondDir(atom1Dir);
        bondDirCounts[firstFromAtom1->getIdx()]+=1;
        atomDirCounts[atom1->getIdx()]+=2;
        atom1ControllingBond=secondFromAtom1;
      }
    } else {
      // dir2 has been set, and dir1 hasn't: we're dealing with a stereochem
      // specification on a ring double bond:
      setFromBond1=false;
      // at least one of the bonds on atom2 has its directionality set already:
      if(bondDirCounts[firstFromAtom2->getIdx()]>0){
        // The second bond's direction has been set at some earlier point:
        atom2Dir = firstFromAtom2->getBondDir();
        bondDirCounts[firstFromAtom2->getIdx()] += 1;
        atomDirCounts[atom2->getIdx()]+=1;
        if(secondFromAtom2){
          // both bonds have their directionalities set, make sure
          // they are compatible:
          if( firstFromAtom2->getBondDir()==secondFromAtom2->getBondDir() &&
              bondDirCounts[firstFromAtom1->getIdx()] ){
            CHECK_INVARIANT(((firstFromAtom2->getBeginAtomIdx()==atom2->getIdx()) ^
                             (secondFromAtom2->getBeginAtomIdx()==atom2->getIdx()))
                            ,"inconsistent state");
          }
        }
      } else {
        // the second bond must be present and setting the direction:
        CHECK_INVARIANT(secondFromAtom2,"inconsistent state");
        CHECK_INVARIANT(bondDirCounts[secondFromAtom2->getIdx()]>0,"inconsistent state");
        // It must be the second bond setting the direction.
        // This happens when the bond dir is set in a branch:
        //        v- this double bond
        //   CC(/C=P/N)=N/O
        //      ^- the second bond sets the direction
        // or when the first bond is a ring closure from an
        // earlier traversed atom:
        //             v- this double bond
        //   NC1=NOC/C1=N\O
        //     ^- this closure ends up being the first bond,
        //        and it does not set the direction.
        //
        // This addresses parts of Issue 185 and sf.net Issue 1842174
        // 
        atom2Dir = secondFromAtom2->getBondDir();

        firstFromAtom2->setBondDir(atom2Dir);
        bondDirCounts[firstFromAtom2->getIdx()]+=1;
        atomDirCounts[atom2->getIdx()]+=2;
        atom2ControllingBond=secondFromAtom2;
      }
      //CHECK_INVARIANT(0,"ring stereochemistry not handled");
    } // end of the ring stereochemistry if

    // now set the directionality on the other side:
    if(setFromBond1){
      if( dblBond->getStereo() == Bond::STEREOE ){
        atom2Dir = atom1Dir;
      } else if( dblBond->getStereo() == Bond::STEREOZ ){
        atom2Dir = (atom1Dir==Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
      }
      CHECK_INVARIANT(atom2Dir != Bond::NONE,"stereo not set");

      // If we're not looking at the bonds used to determine the
      // stereochemistry, we need to flip the setting on the other bond:
      const INT_VECT &stereoAtoms=dblBond->getStereoAtoms();
      if(atom1->getDegree()==3 &&
         std::find(stereoAtoms.begin(),stereoAtoms.end(),
                   static_cast<int>(atom1ControllingBond->getOtherAtomIdx(atom1->getIdx())))
         == stereoAtoms.end() ){
        atom2Dir = (atom2Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
      }
      //std::cerr<<" 0 set bond 2: "<<firstFromAtom2->getIdx()<<" "<<atom2Dir<<std::endl;
      if(atom2->getDegree()==3 &&
         std::find(stereoAtoms.begin(),stereoAtoms.end(),
                   static_cast<int>(firstFromAtom2->getOtherAtomIdx(atom2->getIdx())))
         == stereoAtoms.end() ){
        atom2Dir = (atom2Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
      }
      //std::cerr<<" 1 set bond 2: "<<firstFromAtom2->getIdx()<<" "<<atom2Dir<<std::endl;
      firstFromAtom2->setBondDir(atom2Dir);

      // this block of code is no longer needed
      // if(firstFromAtom2->hasProp(common_properties::_TraversalRingClosureBond)){
      //   // another nice one: we're traversing and come to a ring
      //   // closure bond that has directionality set. This is going to
      //   // have its direction swapped on writing so we need to
      //   // pre-emptively swap it here.
      //   // example situation for this is a non-canonical traversal of
      //   //   C1CCCCN/C=C/1
      //   // starting at atom 0, we hit it on encountering the final bond.
      //   switchBondDir(firstFromAtom2);
      // }

      bondDirCounts[firstFromAtom2->getIdx()] += 1;
      atomDirCounts[atom2->getIdx()]+=1;
    } else {

      // we come before a ring closure:
      if( dblBond->getStereo() == Bond::STEREOZ ){
        atom1Dir = atom2Dir;
      } else if( dblBond->getStereo() == Bond::STEREOE ){
        atom1Dir = (atom2Dir==Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
      }
      CHECK_INVARIANT(atom1Dir != Bond::NONE,"stereo not set");
      // If we're not looking at the bonds used to determine the
      // stereochemistry, we need to flip the setting on the other bond:
      const INT_VECT &stereoAtoms=dblBond->getStereoAtoms();
      if(atom2->getDegree()==3 &&
         std::find(stereoAtoms.begin(),stereoAtoms.end(),
                   static_cast<int>(atom2ControllingBond->getOtherAtomIdx(atom2->getIdx())))
         == stereoAtoms.end() ){
        //std::cerr<<"flip 1"<<std::endl;
        atom1Dir = (atom1Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
      }
      if(atom1->getDegree()==3 &&
         std::find(stereoAtoms.begin(),stereoAtoms.end(),
                   static_cast<int>(firstFromAtom1->getOtherAtomIdx(atom1->getIdx())))
         == stereoAtoms.end() ){
        //std::cerr<<"flip 2"<<std::endl;
        atom1Dir = (atom1Dir == Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
      }
        
      firstFromAtom1->setBondDir(atom1Dir);
      switchBondDir(firstFromAtom1);
      bondDirCounts[firstFromAtom1->getIdx()] += 1;
      atomDirCounts[atom1->getIdx()]+=1;

    }

    // -----------------------------------
    //
    // Check if there are other bonds from atoms 1 and 2 that need
    // to have their directionalities set:
    ///
    if(atom1->getDegree() == 3){
      if(!bondDirCounts[secondFromAtom1->getIdx()] ){
        // This bond (the second bond from the starting atom of the double bond)
        // is a special case.  It's going to appear in a branch in the smiles:
        //     X\C(\Y)=C/Z
        //         ^
        //         |- here
        // so it actually needs to go down with the *same* direction as the
        // bond that's already been set (because "pulling the bond out of the
        // branch" reverses its direction).
        // A quick example.  This SMILES:
        //     F/C(\Cl)=C/F
        // is *wrong*. This is the correct form:
        //     F/C(/Cl)=C/F
        // So, since we want this bond to have the opposite direction to the
        // other one, we put it in with the same direction.
        // This was Issue 183
        secondFromAtom1->setBondDir(firstFromAtom1->getBondDir());
      }
      bondDirCounts[secondFromAtom1->getIdx()] += 1;
      atomDirCounts[atom1->getIdx()]+=1;
    }

    if(atom2->getDegree() == 3){
      if(!bondDirCounts[secondFromAtom2->getIdx()] ){
        // Here we set the bond direction to be opposite the other one (since
        // both come after the atom connected to the double bond).
        Bond::BondDir otherDir;
        if(!secondFromAtom2->hasProp(common_properties::_TraversalRingClosureBond)){
          otherDir = (firstFromAtom2->getBondDir()==Bond::ENDUPRIGHT) ? Bond::ENDDOWNRIGHT : Bond::ENDUPRIGHT;
        } else {
          // another one those irritating little reversal things due to
          // ring closures
          otherDir = firstFromAtom2->getBondDir();
        }
        secondFromAtom2->setBondDir(otherDir);
      }
      bondDirCounts[secondFromAtom2->getIdx()] += 1;
      atomDirCounts[atom2->getIdx()]+=1;
      //std::cerr<<"   other: "<<secondFromAtom2->getIdx()<<" "<<otherDir<<std::endl;
    }

    if(setFromBond1){
      // This is an odd case... The bonds off the beginning atom are
      // after the start atom in the traversal stack.  These need to
      // have their directions reversed.  An example SMILES (unlikely
      // to actually traverse this way is:
      //   C(=C/O)/F    or C(/F)=C/O
      // That bond is Z, without the reversal, this would come out:
      //   C(=C/O)\F    or C(\F)=C/O
      // which is E.
      //
      // In the case of three-coordinate atoms, we don't need to flip
      // the second bond because the Issue 183 fix (above) already got
      // that one.
      //
      // This was Issue 191 and continued into sf.net issue 1842174
      if( bondVisitOrders[atom1ControllingBond->getIdx()] >
          atomVisitOrders[atom1->getIdx()]){
        if(bondDirCounts[atom1ControllingBond->getIdx()]==1){
          if(!atom1ControllingBond->hasProp(common_properties::_TraversalRingClosureBond) ){
            //std::cerr<<"  switcheroo 1"<<std::endl;
            switchBondDir(atom1ControllingBond);
          }
        } else if(bondDirCounts[firstFromAtom2->getIdx()]==1){
          // the controlling bond at atom1 is being set by someone else, flip the direction
          // on the atom2 bond instead:
          //std::cerr<<"  switcheroo 2"<<std::endl;
          switchBondDir(firstFromAtom2);
          if(secondFromAtom2 && bondDirCounts[secondFromAtom2->getIdx()]>=1){
            switchBondDir(secondFromAtom2);
          }
        }
      }
    }

    // something to watch out for here. For this molecule and traversal order:
    //   0 1 2 3  4 5 6  7 8  <- atom numbers
    //   C/C=C/C(/N=C/C)=C/C
    //        ^  ^
    //        |--|-- these two bonds must match in direction or the SMILES
    //               is inconsistent (according to Daylight, Marvin does ok with it)
    // That means that the direction of the bond from atom 3->4 needs to be set
    // when the bond from 2->3 is set.
    // 
    // I believe we only need to worry about this for the bonds from atom2.
    const Atom *atom3=firstFromAtom2->getOtherAtom(atom2);
    if(atom3->getDegree()==3){
      Bond *otherAtom3Bond=NULL;
      bool dblBondPresent=false;
      atomBonds = mol.getAtomBonds(atom3);
      while( atomBonds.first != atomBonds.second ){
        Bond *tbond=mol[*atomBonds.first].get();
        if(tbond->getBondType()==Bond::DOUBLE && tbond->getStereo()>Bond::STEREOANY){
          dblBondPresent=true;
        } else if((tbond->getBondType()==Bond::SINGLE) && (tbond!=firstFromAtom2)) {
          otherAtom3Bond=tbond;
        }
        atomBonds.first++;
      }
      if(dblBondPresent && otherAtom3Bond){
        //std::cerr<<"set!"<<std::endl;
        otherAtom3Bond->setBondDir(firstFromAtom2->getBondDir());
        bondDirCounts[otherAtom3Bond->getIdx()]+=1;
        atomDirCounts[atom3->getIdx()]+=1;
      }
    }
  }

  // finds cycles
  void dfsFindCycles(ROMol &mol,int atomIdx,int inBondIdx,
                     std::vector<AtomColors> &colors,
                     const UINT_VECT &ranks,
                     INT_VECT &atomOrders,
                     VECT_INT_VECT &atomRingClosures,
                     const boost::dynamic_bitset<> *bondsInPlay,
                     const std::vector<std::string> *bondSymbols
                     ){
    Atom *atom = mol.getAtomWithIdx(atomIdx);
    atomOrders.push_back(atomIdx);

    colors[atomIdx] = GREY_NODE;

    // ---------------------
    //
    //  Build the list of possible destinations from here
    //
    // ---------------------
    std::vector< PossibleType > possibles;
    possibles.resize(0);
    ROMol::OBOND_ITER_PAIR bondsPair = mol.getAtomBonds(atom);
    possibles.reserve(bondsPair.second-bondsPair.first);

    while(bondsPair.first != bondsPair.second){
      BOND_SPTR theBond = mol[*(bondsPair.first)];
      bondsPair.first++;
      if(bondsInPlay && !(*bondsInPlay)[theBond->getIdx()]) continue;
      if(inBondIdx<0 || theBond->getIdx() != static_cast<unsigned int>(inBondIdx)){
        int otherIdx = theBond->getOtherAtomIdx(atomIdx);
        long rank=ranks[otherIdx];
        // ---------------------
        //
        // things are a bit more complicated if we are sitting on a
        // ring atom. we would like to traverse first to the
        // ring-closure atoms, then to atoms outside the ring first,
        // then to atoms in the ring that haven't already been visited
        // (non-ring-closure atoms).
        // 
        //  Here's how the black magic works:
        //   - non-ring atom neighbors have their original ranks
        //   - ring atom neighbors have this added to their ranks:
        //       (MAX_BONDTYPE - bondOrder)*MAX_NATOMS*MAX_NATOMS
        //   - ring-closure neighbors lose a factor of:
        //       (MAX_BONDTYPE+1)*MAX_NATOMS*MAX_NATOMS
        //
        //  This tactic biases us to traverse to non-ring neighbors first,
        //  original ordering if bond orders are all equal... crafty, neh?
        //  
        // ---------------------
        if( colors[otherIdx] == GREY_NODE ) {
          rank -= static_cast<int>(MAX_BONDTYPE+1) *
            MAX_NATOMS*MAX_NATOMS;
          if(!bondSymbols){
            rank += static_cast<int>(MAX_BONDTYPE - theBond->getBondType()) *
              MAX_NATOMS;
          } else {
            const std::string &symb=(*bondSymbols)[theBond->getIdx()];
            boost::uint32_t hsh=gboost::hash_range(symb.begin(),symb.end());
            rank += (hsh%MAX_NATOMS) *  MAX_NATOMS;
          }
        } else if( theBond->getOwningMol().getRingInfo()->numBondRings(theBond->getIdx()) ){
          if(!bondSymbols){
            rank += static_cast<int>(MAX_BONDTYPE - theBond->getBondType()) *
              MAX_NATOMS*MAX_NATOMS;
          } else {
            const std::string &symb=(*bondSymbols)[theBond->getIdx()];
            boost::uint32_t hsh=gboost::hash_range(symb.begin(),symb.end());
            rank += (hsh%MAX_NATOMS)*MAX_NATOMS*MAX_NATOMS;
          }
        }
       //std::cerr<<"aIdx: "<< atomIdx <<"   p: "<<otherIdx<<" Rank: "<<ranks[otherIdx] <<" "<<colors[otherIdx]<<" "<<theBond->getBondType()<<" "<<rank<<std::endl;
        possibles.push_back(PossibleType(rank,otherIdx,theBond.get()));
      }
    }

    // ---------------------
    //
    //  Sort on ranks
    //
    // ---------------------
    std::sort(possibles.begin(),possibles.end(),_possibleCompare());


    // ---------------------
    //
    //  Now work the children
    //
    // ---------------------
    for(std::vector<PossibleType>::iterator possiblesIt=possibles.begin();
        possiblesIt!=possibles.end();
        possiblesIt++){
      int possibleIdx = possiblesIt->get<1>();
      Bond *bond = possiblesIt->get<2>();
      Atom *otherAtom=mol.getAtomWithIdx(possibleIdx);
      switch(colors[possibleIdx]){
      case WHITE_NODE:
        // -----
        // we haven't seen this node at all before, traverse
        // -----
        dfsFindCycles(mol,possibleIdx,bond->getIdx(),colors,
                    ranks,atomOrders,
                    atomRingClosures,
                    bondsInPlay,bondSymbols);
        break;
      case GREY_NODE:
        // -----
        // we've seen this, but haven't finished it (we're finishing a ring)
        // -----
        atomRingClosures[possibleIdx].push_back(bond->getIdx());
        atomRingClosures[atomIdx].push_back(bond->getIdx());
        break;
      default:
        // -----
        // this node has been finished. don't do anything.
        // -----
        break;
      }
    }
    colors[atomIdx] = BLACK_NODE;
  }

  void dfsBuildStack(ROMol &mol,int atomIdx,int inBondIdx,
                     std::vector<AtomColors> &colors,
                     VECT_INT_VECT &cycles,
                     const UINT_VECT &ranks,
                     INT_VECT &cyclesAvailable,
                     MolStack &molStack,
                     INT_VECT &atomOrders,
                     INT_VECT &bondVisitOrders,
                     VECT_INT_VECT &atomRingClosures,
                     std::vector<INT_LIST> &atomTraversalBondOrder,
                     const boost::dynamic_bitset<> *bondsInPlay,
                     const std::vector<std::string> *bondSymbols
                     ){
#if 0
    std::cerr<<"traverse from atom: "<<atomIdx<<" via bond "<<inBondIdx<<" num cycles available: "
             <<std::count(cyclesAvailable.begin(),cyclesAvailable.end(),1)<<std::endl;
#endif
    Atom *atom = mol.getAtomWithIdx(atomIdx);
    INT_LIST directTravList,cycleEndList;
    boost::dynamic_bitset<> seenFromHere(mol.getNumAtoms());
    
    seenFromHere.set(atomIdx);
    molStack.push_back(MolStackElem(atom));
    atomOrders[atom->getIdx()] = molStack.size();
    colors[atomIdx] = GREY_NODE;

    INT_LIST travList;
    if(inBondIdx>=0) travList.push_back(inBondIdx);

    
    // ---------------------
    //
    //  Add any ring closures
    //
    // ---------------------
    if(atomRingClosures[atomIdx].size()){
      std::vector<unsigned int> ringsClosed;
      BOOST_FOREACH(int bIdx,atomRingClosures[atomIdx]){
        travList.push_back(bIdx);
        Bond *bond = mol.getBondWithIdx(bIdx);
        seenFromHere.set(bond->getOtherAtomIdx(atomIdx));
        unsigned int ringIdx;
        if(bond->getPropIfPresent(common_properties::_TraversalRingClosureBond, ringIdx)){
          // this is end of the ring closure
          // we can just pull the ring index from the bond itself:
          molStack.push_back(MolStackElem(bond,atomIdx));
          bondVisitOrders[bIdx]=molStack.size();
          molStack.push_back(MolStackElem(ringIdx));
          // don't make the ring digit immediately available again: we don't want to have the same
          // ring digit opening and closing rings on an atom.
          ringsClosed.push_back(ringIdx-1);
        } else {
          // this is the beginning of the ring closure, we need to come up with a ring index:
          INT_VECT::const_iterator cAIt=std::find(cyclesAvailable.begin(),
                                                  cyclesAvailable.end(),1);
          if(cAIt==cyclesAvailable.end()){
            throw ValueErrorException("Too many rings open at once. SMILES cannot be generated.");
          }
          unsigned int lowestRingIdx =  cAIt-cyclesAvailable.begin();
          cyclesAvailable[lowestRingIdx] = 0;
          ++lowestRingIdx;
          bond->setProp(common_properties::_TraversalRingClosureBond,lowestRingIdx);
          molStack.push_back(MolStackElem(lowestRingIdx));
        }
      }
      BOOST_FOREACH(unsigned int ringIdx,ringsClosed){
        cyclesAvailable[ringIdx]=1;
      }
    }

    
    // ---------------------
    //
    //  Build the list of possible destinations from here
    //
    // ---------------------
    std::vector< PossibleType > possibles;
    possibles.resize(0);
    ROMol::OBOND_ITER_PAIR bondsPair = mol.getAtomBonds(atom);
    possibles.reserve(bondsPair.second-bondsPair.first);

    while(bondsPair.first != bondsPair.second){
      BOND_SPTR theBond = mol[*(bondsPair.first)];
      bondsPair.first++;
      if(bondsInPlay && !(*bondsInPlay)[theBond->getIdx()]) continue;
      if(inBondIdx<0 || theBond->getIdx() != static_cast<unsigned int>(inBondIdx)){
        int otherIdx = theBond->getOtherAtomIdx(atomIdx);
        // ---------------------
        //
        // This time we skip the ring-closure atoms (we did them
        // above); we want to traverse first to atoms outside the ring
        // then to atoms in the ring that haven't already been visited
        // (non-ring-closure atoms).
        // 
        // otherwise it's the same ranking logic as above
        //  
        // ---------------------
        if( colors[otherIdx] != WHITE_NODE || seenFromHere[otherIdx] ) {
          // ring closure or finished atom... skip it.
          continue;
        }
        unsigned long rank=ranks[otherIdx];
        if( theBond->getOwningMol().getRingInfo()->numBondRings(theBond->getIdx()) ){
          if(!bondSymbols){
            rank += static_cast<int>(MAX_BONDTYPE - theBond->getBondType()) *
              MAX_NATOMS*MAX_NATOMS;
          } else {
            const std::string &symb=(*bondSymbols)[theBond->getIdx()];
            boost::uint32_t hsh=gboost::hash_range(symb.begin(),symb.end());
            rank += (hsh%MAX_NATOMS)*MAX_NATOMS*MAX_NATOMS;
          }
        }
        //std::cerr<<"   p: "<<otherIdx<<" "<<colors[otherIdx]<<" "<<theBond->getBondType()<<" "<<rank<<std::endl;
        possibles.push_back(PossibleType(rank,otherIdx,theBond.get()));
      }
    }

    // ---------------------
    //
    //  Sort on ranks
    //
    // ---------------------
    std::sort(possibles.begin(),possibles.end(),_possibleCompare());

    // ---------------------
    //
    //  Now work the children
    //
    // ---------------------
    for(std::vector<PossibleType>::iterator possiblesIt=possibles.begin();
        possiblesIt!=possibles.end();
        possiblesIt++){
      int possibleIdx = possiblesIt->get<1>();
      if(colors[possibleIdx]!=WHITE_NODE){
        // we're either done or it's a ring-closure, which we already processed...
        // this test isn't strictly required, because we only added WHITE notes to
        // the possibles list, but it seems logical to document it
        continue;
      }
      Bond *bond = possiblesIt->get<2>();
      Atom *otherAtom=mol.getAtomWithIdx(possibleIdx);
      unsigned int lowestRingIdx;
      INT_VECT::const_iterator cAIt;
      // ww might have some residual data from earlier calls, clean that up:
      if(otherAtom->hasProp(common_properties::_TraversalBondIndexOrder)){
        otherAtom->clearProp(common_properties::_TraversalBondIndexOrder);
      }
      travList.push_back(bond->getIdx());
      if(possiblesIt+1 != possibles.end()){
        // we're branching
        molStack.push_back(MolStackElem("(",possiblesIt-possibles.begin()));
      }
      molStack.push_back(MolStackElem(bond,atomIdx));
      bondVisitOrders[bond->getIdx()]=molStack.size();
      dfsBuildStack(mol,possibleIdx,bond->getIdx(),colors,
                    cycles,ranks,cyclesAvailable,molStack,
                    atomOrders,bondVisitOrders,atomRingClosures,atomTraversalBondOrder,
                    bondsInPlay,bondSymbols);
      if(possiblesIt+1 != possibles.end()){
        molStack.push_back(MolStackElem(")",possiblesIt-possibles.begin()));
      }
    }

    atomTraversalBondOrder[atom->getIdx()]=travList;
    colors[atomIdx] = BLACK_NODE;
  }

  
  void canonicalDFSTraversal(ROMol &mol,int atomIdx,int inBondIdx,
                             std::vector<AtomColors> &colors,
                             VECT_INT_VECT &cycles,
                             const UINT_VECT &ranks,
                             INT_VECT &cyclesAvailable,
                             MolStack &molStack,
                             INT_VECT &atomOrders,
                             INT_VECT &bondVisitOrders,
                             VECT_INT_VECT &atomRingClosures,
                             std::vector<INT_LIST> &atomTraversalBondOrder,
                             const boost::dynamic_bitset<> *bondsInPlay,
                             const std::vector<std::string> *bondSymbols
                             ){
    PRECONDITION(colors.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(ranks.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(atomOrders.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(bondVisitOrders.size()>=mol.getNumBonds(),"vector too small");
    PRECONDITION(atomRingClosures.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(atomTraversalBondOrder.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(!bondsInPlay || bondsInPlay->size()>=mol.getNumBonds(),"bondsInPlay too small");
    PRECONDITION(!bondSymbols || bondSymbols->size()>=mol.getNumBonds(),"bondSymbols too small");

    std::vector<AtomColors> tcolors;
    tcolors.resize(colors.size());
    std::copy(colors.begin(),colors.end(),tcolors.begin());
    dfsFindCycles(mol,atomIdx,inBondIdx,
                  tcolors,
                  ranks,
                  atomOrders,
                  atomRingClosures,
                  bondsInPlay,
                  bondSymbols);
    dfsBuildStack(mol,atomIdx,inBondIdx,
                  colors,
                  cycles,
                  ranks,
                  cyclesAvailable,
                  molStack,
                  atomOrders,
                  bondVisitOrders,
                  atomRingClosures,
                  atomTraversalBondOrder,
                  bondsInPlay,
                  bondSymbols);

  }

  bool canHaveDirection(const Bond *bond){
    PRECONDITION(bond,"bad bond");
    Bond::BondType bondType= bond->getBondType();
    return (bondType==Bond::SINGLE || bondType==Bond::AROMATIC);
  }

  void clearBondDirs(ROMol &mol,Bond *refBond,const Atom *fromAtom,
                     INT_VECT &bondDirCounts,
                     INT_VECT &atomDirCounts,
                     const INT_VECT &bondVisitOrders){
    PRECONDITION(bondDirCounts.size()>=mol.getNumBonds(),"bad dirCount size");
    PRECONDITION(refBond,"bad bond");
    PRECONDITION(&refBond->getOwningMol()==&mol,"bad bond");
    PRECONDITION(fromAtom,"bad atom");
    PRECONDITION(&fromAtom->getOwningMol()==&mol,"bad bond");

#if 0
    std::copy(bondDirCounts.begin(),bondDirCounts.end(),std::ostream_iterator<int>(std::cerr,", "));
    std::cerr<<"\n";
    std::copy(atomDirCounts.begin(),atomDirCounts.end(),std::ostream_iterator<int>(std::cerr,", "));
    std::cerr<<"\n";
    std::cerr<<"cBD: bond: "<<refBond->getIdx()<<" atom: "<<fromAtom->getIdx()<<": ";
#endif
    ROMol::OEDGE_ITER beg,end;
    boost::tie(beg,end) = mol.getAtomBonds(fromAtom);
    bool nbrPossible=false,adjusted=false;
    while(beg!=end){
      Bond *oBond=mol[*beg].get();
      //std::cerr<<"  >>"<<oBond->getIdx()<<" "<<canHaveDirection(oBond)<<" "<<bondDirCounts[oBond->getIdx()]<<"-"<<bondDirCounts[refBond->getIdx()]<<" "<<atomDirCounts[oBond->getBeginAtomIdx()]<<"-"<<atomDirCounts[oBond->getEndAtomIdx()]<<std::endl;
      if( oBond != refBond && canHaveDirection(oBond) ){
        nbrPossible=true;
        if((bondDirCounts[oBond->getIdx()] >= bondDirCounts[refBond->getIdx()])
           &&
           atomDirCounts[oBond->getBeginAtomIdx()]!=1 && atomDirCounts[oBond->getEndAtomIdx()]!=1){
          adjusted=true;
          bondDirCounts[oBond->getIdx()] -= 1;
          if(!bondDirCounts[oBond->getIdx()]){
            // no one is setting the direction here:
            oBond->setBondDir(Bond::NONE);
            atomDirCounts[oBond->getBeginAtomIdx()]-=1;
            atomDirCounts[oBond->getEndAtomIdx()]-=1;
            //std::cerr<<"ob:"<<oBond->getIdx()<<" ";
          }
        }
      }
      beg++;
    }
    if(nbrPossible && !adjusted &&
       atomDirCounts[refBond->getBeginAtomIdx()]!=1 && atomDirCounts[refBond->getEndAtomIdx()]!=1){
      // we found a neighbor that could have directionality set,
      // but it had a lower bondDirCount than us, so we must
      // need to be adjusted:
      bondDirCounts[refBond->getIdx()] -= 1;
      if(!bondDirCounts[refBond->getIdx()]){
        refBond->setBondDir(Bond::NONE);
        atomDirCounts[refBond->getBeginAtomIdx()]-=1;
        atomDirCounts[refBond->getEndAtomIdx()]-=1;
        //std::cerr<<"rb:"<<refBond->getIdx()<<" ";
      }
    }
    //std::cerr<<std::endl;
  }

  void removeRedundantBondDirSpecs(ROMol &mol,MolStack &molStack,INT_VECT &bondDirCounts,INT_VECT &atomDirCounts,
                                   const INT_VECT &bondVisitOrders){
    PRECONDITION(bondDirCounts.size()>=mol.getNumBonds(),"bad dirCount size");
#if 0
    std::cerr<<"rRBDS: ";
    mol.debugMol(std::cerr);
    std::copy(bondDirCounts.begin(),bondDirCounts.end(),std::ostream_iterator<int>(std::cerr,", "));
    std::cerr<<"\n";
#endif    
    // find bonds that have directions indicated that are redundant:
    for(MolStack::iterator msI=molStack.begin();
        msI!=molStack.end(); msI++) {
      if( msI->type == MOL_STACK_BOND ){
        Bond *tBond = msI->obj.bond;
        const Atom *canonBeginAtom=mol.getAtomWithIdx(msI->number);
        const Atom *canonEndAtom=mol.getAtomWithIdx(tBond->getOtherAtomIdx(msI->number));
        if(canHaveDirection(tBond) &&
           bondDirCounts[tBond->getIdx()]>=1 ) {
          // start by finding the double bond that sets tBond's direction:
          const Atom *dblBondAtom=NULL;
          ROMol::OEDGE_ITER beg,end;
          boost::tie(beg,end) = mol.getAtomBonds(canonBeginAtom);
          while(beg!=end){
            if( mol[*beg].get() != tBond && mol[*beg]->getBondType()==Bond::DOUBLE &&
                mol[*beg]->getStereo() > Bond::STEREOANY ){
              dblBondAtom = canonBeginAtom;//tBond->getOtherAtom(canonBeginAtom);
              break;
            }
            beg++;
          }
          if(dblBondAtom != NULL){
            clearBondDirs(mol,tBond,dblBondAtom,bondDirCounts,atomDirCounts,bondVisitOrders);
          }
          dblBondAtom = NULL;
          boost::tie(beg,end) = mol.getAtomBonds(canonEndAtom);
          while(beg!=end){
            if( mol[*beg].get() != tBond && mol[*beg]->getBondType()==Bond::DOUBLE &&
                mol[*beg]->getStereo() > Bond::STEREOANY ){
              dblBondAtom = canonEndAtom;//tBond->getOtherAtom(canonEndAtom);
              break;
            }
            beg++;
          }
          if(dblBondAtom != NULL){
            clearBondDirs(mol,tBond,dblBondAtom,bondDirCounts,atomDirCounts,bondVisitOrders);
          }
        } else if(tBond->getBondDir()!=Bond::NONE){
          // we aren't supposed to have a direction set, but we do:
          tBond->setBondDir(Bond::NONE);
        }
      }
    }
  }

  void canonicalizeFragment(ROMol &mol,int atomIdx,
                            std::vector<AtomColors> &colors,
                            const UINT_VECT &ranks,
                            MolStack &molStack,
                            const boost::dynamic_bitset<> *bondsInPlay,
                            const std::vector<std::string> *bondSymbols,
                            bool doIsomericSmiles){
    PRECONDITION(colors.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(ranks.size()>=mol.getNumAtoms(),"vector too small");
    PRECONDITION(!bondsInPlay || bondsInPlay->size()>=mol.getNumBonds(),"bondsInPlay too small");
    PRECONDITION(!bondSymbols || bondSymbols->size()>=mol.getNumBonds(),"bondSymbols too small");
    int nAtoms=mol.getNumAtoms();
    
    INT_VECT atomVisitOrders(nAtoms,0);
    INT_VECT bondVisitOrders(mol.getNumBonds(),0);
    INT_VECT bondDirCounts(mol.getNumBonds(),0);
    INT_VECT atomDirCounts(nAtoms,0);
    INT_VECT cyclesAvailable(MAX_CYCLES,1);

    VECT_INT_VECT cycles(nAtoms);
    for(VECT_INT_VECT_I vviIt=cycles.begin();vviIt!=cycles.end();++vviIt) vviIt->resize(0);

    boost::dynamic_bitset<> ringStereoChemAdjusted(nAtoms);
    
    // make sure that we've done the stereo perception:
    if(!mol.hasProp(common_properties::_StereochemDone)){
      MolOps::assignStereochemistry(mol,false);
    }

    // we need ring information; make sure findSSSR has been called before
    // if not call now
    if ( !mol.getRingInfo()->isInitialized() ) {
      MolOps::findSSSR(mol);
    }
    mol.getAtomWithIdx(atomIdx)->setProp("_TraversalStartPoint",true);

    VECT_INT_VECT atomRingClosures(nAtoms);
    std::vector<INT_LIST> atomTraversalBondOrder(nAtoms);
    Canon::canonicalDFSTraversal(mol,atomIdx,-1,colors,cycles,
                                 ranks,cyclesAvailable,molStack,atomVisitOrders,
                                 bondVisitOrders,atomRingClosures,atomTraversalBondOrder,
                                 bondsInPlay,bondSymbols);

    PRECONDITION(!molStack.empty(), "Empty stack.");
    PRECONDITION(molStack.begin()->type == MOL_STACK_ATOM, "Corrupted stack. First element should be an atom.");

    // collect some information about traversal order on chiral atoms
    bool *numSwapsChiralAtoms=(bool *)malloc(nAtoms*sizeof(bool));
    memset(numSwapsChiralAtoms,0,nAtoms*sizeof(bool));
    if(doIsomericSmiles){
      for(ROMol::AtomIterator atomIt=mol.beginAtoms();atomIt!=mol.endAtoms();++atomIt){
        if((*atomIt)->getChiralTag()!=Atom::CHI_UNSPECIFIED){
          ROMol::OEDGE_ITER beg,end;
          boost::tie(beg,end) = mol.getAtomBonds(*atomIt);
          while(beg!=end){
            if(bondsInPlay && !(*bondsInPlay)[mol[*beg]->getIdx()]){
              (*atomIt)->setProp(common_properties::_brokenChirality,true);
              break;
            }
            ++beg;
          }
          if((*atomIt)->hasProp(common_properties::_brokenChirality)){
            continue;
          }
          INT_LIST trueOrder = atomTraversalBondOrder[(*atomIt)->getIdx()];
          //Test if the atom is in current fragment
          if(trueOrder.size()>0){
            int nSwaps=  (*atomIt)->getPerturbationOrder(trueOrder);
            if((*atomIt)->getDegree()==3 && molStack.begin()->obj.atom->getIdx() == (*atomIt)->getIdx()){
              // This is a special case. Here's an example:
              //   Our internal representation of a chiral center is equivalent to:
              //     [C@](F)(O)(C)[H]
              //   we'll be dumping it without the H, which entails a reordering:
              //     [C@@H](F)(O)C
              ++nSwaps;
            }
            if(nSwaps%2){
              numSwapsChiralAtoms[(*atomIt)->getIdx()] = 1;
            }
          }
        }
      }
    }

    // remove the current directions on single bonds around double bonds:
    for(ROMol::BondIterator bondIt=mol.beginBonds();
        bondIt!=mol.endBonds();
        ++bondIt){
      Bond::BondDir dir = (*bondIt)->getBondDir();
      if (dir == Bond::ENDDOWNRIGHT || dir == Bond::ENDUPRIGHT) {
        (*bondIt)->setBondDir(Bond::NONE);
      }
    }

#if 0
    std::cerr<<"<11111111"<<std::endl;

    std::cerr<<"----------------------------------------->"<<std::endl;
    mol.debugMol(std::cerr);
#endif


    //std::cerr<<"----->\ntraversal stack:"<<std::endl;
    // traverse the stack and canonicalize double bonds and atoms with (ring) stereochemistry
    for(MolStack::iterator msI=molStack.begin();
        msI!=molStack.end(); ++msI){
#if 0
      if(msI->type == MOL_STACK_ATOM) std::cerr<<" atom: "<<msI->obj.atom->getIdx()<<std::endl;
      else if(msI->type == MOL_STACK_BOND) std::cerr<<" bond: "<<msI->obj.bond->getIdx()<<" "<<msI->number<<" "<<msI->obj.bond->getBeginAtomIdx()<<"-"<<msI->obj.bond->getEndAtomIdx()<<" order: "<<msI->obj.bond->getBondType()<<std::endl;
      else if(msI->type == MOL_STACK_RING) std::cerr<<" ring: "<<msI->number<<std::endl;
      else if(msI->type == MOL_STACK_BRANCH_OPEN) std::cerr<<" branch open"<<std::endl;
      else if(msI->type == MOL_STACK_BRANCH_CLOSE) std::cerr<<" branch close"<<std::endl;
#endif      
      if(msI->type == MOL_STACK_BOND &&
          msI->obj.bond->getBondType() == Bond::DOUBLE &&
          msI->obj.bond->getStereo() > Bond::STEREOANY){
        if(msI->obj.bond->getStereoAtoms().size()>=2){
          Canon::canonicalizeDoubleBond(msI->obj.bond,bondVisitOrders,atomVisitOrders,
              bondDirCounts,atomDirCounts);
        } else {
          // bad stereo spec:
          msI->obj.bond->setStereo(Bond::STEREONONE);
        }
      }
      if(doIsomericSmiles){
        if(msI->type == MOL_STACK_ATOM && msI->obj.atom->getChiralTag()!=Atom::CHI_UNSPECIFIED
            && !msI->obj.atom->hasProp(common_properties::_brokenChirality)){
          if(msI->obj.atom->hasProp(common_properties::_ringStereoAtoms)){
            if(!ringStereoChemAdjusted[msI->obj.atom->getIdx()]){
              msI->obj.atom->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
              ringStereoChemAdjusted.set(msI->obj.atom->getIdx());
            }
            const INT_VECT &ringStereoAtoms=msI->obj.atom->getProp<INT_VECT>(common_properties::_ringStereoAtoms);
            BOOST_FOREACH(int nbrV,ringStereoAtoms){
              int nbrIdx=abs(nbrV)-1;
              //Adjust the chiraliy flag of the ring stereo atoms according to the first one
              if(!ringStereoChemAdjusted[nbrIdx] &&
                  atomVisitOrders[nbrIdx]>atomVisitOrders[msI->obj.atom->getIdx()]){
                mol.getAtomWithIdx(nbrIdx)->setChiralTag(msI->obj.atom->getChiralTag());
                if(nbrV<0){
                  mol.getAtomWithIdx(nbrIdx)->invertChirality();
                }
                //Odd number of swaps for first chiral ring atom --> needs to be swapped but we want to retain chirality
                if(numSwapsChiralAtoms[msI->obj.atom->getIdx()]){
                  //Odd number of swaps for chiral ring neighbor --> needs to be swapped but we want to retain chirality
                  if(!numSwapsChiralAtoms[nbrIdx]){
                    mol.getAtomWithIdx(nbrIdx)->invertChirality();
                  }
                }
                //Even number of swaps for first chiral ring atom --> don't need to be swapped
                else{
                  //Odd number of swaps for chiral ring neighbor --> needs to be swapped
                  if(numSwapsChiralAtoms[nbrIdx]){
                    mol.getAtomWithIdx(nbrIdx)->invertChirality();
                  }
                }
                ringStereoChemAdjusted.set(nbrIdx);
              }
            }
          }
          else{
            if(msI->obj.atom->getChiralTag()== Atom::CHI_TETRAHEDRAL_CW){
              if((numSwapsChiralAtoms[msI->obj.atom->getIdx()])){
                msI->obj.atom->invertChirality();
              }
            }
            else if(msI->obj.atom->getChiralTag()== Atom::CHI_TETRAHEDRAL_CCW){
              if((numSwapsChiralAtoms[msI->obj.atom->getIdx()])){
                msI->obj.atom->invertChirality();
              }
            }
          }
        }
      }
    }
#if 0
    std::cerr<<"<-----"<<std::endl;

    std::cerr<<"----------------------------------------->"<<std::endl;
    mol.debugMol(std::cerr);
#endif
    Canon::removeRedundantBondDirSpecs(mol,molStack,bondDirCounts,atomDirCounts,bondVisitOrders);
#if 0
    std::cerr<<"----------------------------------------->"<<std::endl;
    mol.debugMol(std::cerr);
    std::cerr<<"----------------------------------------->"<<std::endl;
#endif
    free(numSwapsChiralAtoms);
  }
}
}