C++ OBBitVec类代码示例

  static int DetermineFRJ(OBMol &mol)
  {
    vector<vector<int> >::iterator i;
    vector<vector<int> > cfl;
    //find all continuous graphs in the mol area
    mol.ContigFragList(cfl);

    if (cfl.empty())
      return(0);
    if (cfl.size() == 1)
      return(mol.NumBonds() - mol.NumAtoms() + 1);

    //count up the atoms and bonds belonging to each graph
    OBBond *bond;
    vector<OBBond*>::iterator j;
    int numatoms,numbonds,frj=0;
    OBBitVec frag;
    for (i = cfl.begin();i != cfl.end();++i)
      {
        frag.Clear();
        frag.FromVecInt(*i);
        numatoms = (*i).size();
        numbonds=0;
        for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))
          if (frag.BitIsOn(bond->GetBeginAtomIdx()) &&
              frag.BitIsOn(bond->GetEndAtomIdx()))
            numbonds++;
        frj += numbonds - numatoms + 1;
      }

    return(frj);
  }

  bool OBRingSearch::SaveUniqueRing(deque<int> &d1,deque<int> &d2)
  {
    vector<int> path;
    OBBitVec bv;
    deque<int>::iterator i;

    for (i = d1.begin();i != d1.end();++i)
      {
        bv.SetBitOn(*i);
        path.push_back(*i);
      }

    for (i = d2.begin();i != d2.end();++i)
      {
        bv.SetBitOn(*i);
        path.push_back(*i);
      }

    vector<OBRing*>::iterator j;
    for (j = _rlist.begin();j != _rlist.end();++j)
      if (bv == (*j)->_pathset)
        return(false);

    OBRing *ring = new OBRing(path, bv);
    _rlist.push_back(ring);

    return(true);
  }

  vector<OBBond*> OBResidue::GetBonds(bool exterior) const
  {
    OBAtom         *atom;
    vector<OBBond*> bonds;
    OBBitVec        idxs;
    unsigned int    sz;

    sz = (unsigned int) _atoms.size();
    for ( unsigned int i = 0 ; i < sz ; ++i )
      {
        atom = _atoms[i];
        OBBond *bond;
        vector<OBBond*>::iterator b;
        for (bond = atom->BeginBond(b) ; bond ; bond = atom->NextBond(b))
          {
            if (!idxs.BitIsOn(bond->GetIdx()))
              {
                if (!exterior)
                  {
                    if (bond->GetNbrAtom(atom)->GetResidue() == this)
                      bonds.push_back(&(*bond));
                  }
                else
                  bonds.push_back(&(*bond));

                idxs.SetBitOn(bond->GetIdx());
              }
          }
      }

    return bonds;
  }

 /** Return true if \p bv1 i less than \p bv2
     Lexicographical order, with bit vectors written LSB first.
     \param[in] bv1 A bit vector
     \param[in] bv2 Another bit vector
     \return true if equal, false otherwise
 */
 bool operator< (const OBBitVec & bv1, const OBBitVec & bv2)
 {
   bool rtn = false;
   int next_bit_1 = bv1.NextBit(-1);
   int next_bit_2 = bv2.NextBit(-1);
   bool should_continue = true;
   while (should_continue)
     {
     should_continue = false;
     if (next_bit_1 == -1)
       rtn = (next_bit_2 == -1 ? false : true);
     else if (next_bit_2 == -1)
       rtn = false;
     else if (next_bit_2 < next_bit_1)
       rtn = true;
     else if (next_bit_1 < next_bit_2)
       rtn = false;
     else
       {
       next_bit_1 = bv1.NextBit(next_bit_1);
       next_bit_2 = bv2.NextBit(next_bit_2);
       should_continue = true;
       }
     }
   return rtn;
 }

  void OBRingSearch::RemoveRedundant(int frj)
  {
    int i,j;

    //remove identical rings
    for (i = _rlist.size()-1;i > 0;i--)
      for (j = i-1;j >= 0;j--)
        if ((_rlist[i])->_pathset == (_rlist[j])->_pathset)
          {
            delete _rlist[i];
            _rlist.erase(_rlist.begin()+i);
            break;
          }

    if (_rlist.size() == 0)
      return; // nothing to do

    // handle LSSR
    if (frj < 0) {
      OBMol *mol = _rlist[0]->GetParent();
      std::vector<OBRing*> rlist, rignored;
      for (unsigned int i = 0; i < _rlist.size(); ++i) {
        visitRing(mol, _rlist[i], rlist, rignored);
      }
      for (unsigned int i = 0; i < rignored.size(); ++i)
        delete rignored[i];
      _rlist = rlist;
      return;
    }

    // exit if we already have frj rings
    if (_rlist.size() == (unsigned)frj)
      return;

    //make sure tmp is the same size as the rings
    OBBitVec tmp;
    for (j = 0;j < (signed)_rlist.size();++j)
      tmp = (_rlist[j])->_pathset;

    //remove larger rings that cover the same atoms as smaller rings
    for (i = _rlist.size()-1;i >= 0;i--)
      {
        tmp.Clear();
        for (j = 0;j < (signed)_rlist.size();++j)
          if ((_rlist[j])->_path.size() <= (_rlist[i])->_path.size() && i != j)
            tmp |= (_rlist[j])->_pathset;

        tmp = tmp & (_rlist[i])->_pathset;

        if (tmp == (_rlist[i])->_pathset)
          {
            delete _rlist[i];
            _rlist.erase(_rlist.begin()+i);
          }

        if (_rlist.size() == (unsigned)frj)
          break;
      }
  }

void testAutomorphismMask() {
  // read file: 3 6-rings
  //
  //     /\ /\ /\
  //    |  |  |  |
  //     \/ \/ \/
  //
  cout <<  "testAutomorphismMask" << endl;
  OBMol mol;
  OBConversion conv;
  conv.SetInFormat("cml");
  std::ifstream ifs(OBTestUtil::GetFilename("isomorphism1.cml").c_str());
  OB_REQUIRE( ifs );
  conv.Read(&mol, &ifs);

  OBIsomorphismMapper::Mappings maps;

  // First of all, how many automorphisms are there without any mask?
  // This takes about 20 seconds, so you may want to comment this out while debugging
  FindAutomorphisms(&mol, maps);
  cout << maps.size() << endl;
  OB_ASSERT( maps.size() == 4 );

  // Now, let's remove the bridge (atomId 6) of the central ring.
  //
  //     /\ /\ /\
  //    |  |  |  |
  //     \/    \/
  // both rings can be flipped around exocyclic bond, the whole molecule can be mirrored
  // horizontally, this results in 2 x 2 x 2 = 8 automorphisms
  OBBitVec mask;
  mask.SetRangeOn(1, mol.NumAtoms());
  mask.SetBitOff(6+1);
  FindAutomorphisms(&mol, maps, mask);
  cout << maps.size() << endl;
  for (unsigned int i = 0; i < maps.size(); ++i) {
    OBIsomorphismMapper::Mapping::const_iterator j;
    for (j = maps[i].begin(); j != maps[i].end(); ++j)
      cout << j->second << " ";
    cout << endl;
  }
  OB_ASSERT( maps.size() == 8 );

  // Verify that atom Id 6 does not occur anywhere in the mappings
  OBIsomorphismMapper::Mappings::const_iterator a;
  OBIsomorphismMapper::Mapping::const_iterator b;
  for (a = maps.begin(); a != maps.end(); ++a)
    for (b = a->begin(); b!= a->end(); ++b) {
      OB_ASSERT( b->first != 6 );
      OB_ASSERT( b->second != 6 );
    }
}

  /**
   * This function finds the LSSR containing all relevant cycles. A cycle is
   * relevant if it belongs to at least one minimum cycle basis. Another
   * description is more useful though:
   *
   * A cycle (C) is relevant if:
   * - no smaller cycles C_i, ..., C_k exist such that C = C_1 + ... + C_k
   * - both bonds & atoms are checked
   *
   * This is based on lemma 1 from:
   *
   * P. Vismara, Union of all the minimum cycle bases of a graph, The electronic
   * journal of combinatorics, Vol. 4, 1997
   * http://www.emis.de/journals/EJC/Volume_4/PostScriptfiles/v4i1r9.ps
   */
  void visitRing(OBMol *mol, OBRing *ring, std::vector<OBRing*> &rlist, std::vector<OBRing*> &rignored)
  {
    const std::vector<int> &atoms = ring->_path;
    OBBitVec mask;
    // Make sure mask is the same size as the maximum ring atom/bond index.
    mask.SetBitOn(mol->NumAtoms());
    mask.SetBitOn(mol->NumBonds());

    //
    // Remove larger rings that cover the same atoms as smaller rings.
    //
    mask.Clear();
    for (unsigned int j = 0; j < rlist.size(); ++j)
      // Here we select only smaller rings.
      if (rlist[j]->_path.size() < ring->_path.size())
        mask |= rlist[j]->_pathset;

    mask = mask & ring->_pathset;

    bool containsSmallerAtomRing = (mask == ring->_pathset) ? true : false;

    // Translate ring atom indexes to ring bond indexes.
    std::vector<unsigned int> bonds = atomRingToBondRing(mol, ring->_path);
    OBBitVec bondset;
    for (unsigned int i = 0; i < bonds.size(); ++i)
      bondset.SetBitOn(bonds[i]);

    //
    // Remove larger rings that cover the same bonds as smaller rings.
    //
    mask.Clear();
    for (unsigned int j = 0; j < rlist.size(); ++j) {
      std::vector<unsigned int> otherBonds = atomRingToBondRing(mol, rlist[j]->_path);
      OBBitVec bs;
      for (unsigned int i = 0; i < otherBonds.size(); ++i)
        bs.SetBitOn(otherBonds[i]);

      // Here we select only smaller rings.
      if (otherBonds.size() < bonds.size())
        mask |= bs;
    }

    mask = mask & bondset;

    bool containsSmallerBondRing = (mask == bondset) ? true : false;

    // The ring is part of the LSSR if all it's atoms and bonds are not
    // found in smaller rings.
    if (!containsSmallerAtomRing || !containsSmallerBondRing) {
      rlist.push_back(ring);
    } else {
      rignored.push_back(ring);
    }
  }

      bool isInTerminalSet(const std::vector<unsigned int> &depths,
          const OBBitVec &path, std::size_t i)
      {
        if (!depths[i])
          return false;

        if (path.BitIsSet(i))
          return false;

        return true;
      }

  static void FindRings(OBMol &mol,vector<int> &path,OBBitVec &avisit,
                        OBBitVec &bvisit, int natom,int depth )
  {
    OBAtom *atom;
    OBBond *bond;
    vector<OBBond*>::iterator k;

    // don't return if all atoms are visited
    // (For example, some atoms are in multiple rings!) -GRH
      
    if (avisit[natom])
      {
        int j = depth-1;
        bond=mol.GetBond(path[j--]);
        bond->SetInRing();
        while( j >= 0 )
          {
            bond=mol.GetBond(path[j--]);
            bond->SetInRing();
            (bond->GetBeginAtom())->SetInRing();
            (bond->GetEndAtom())->SetInRing();
            if(bond->GetBeginAtomIdx()==static_cast<unsigned int>(natom) || bond->
               GetEndAtomIdx()==static_cast<unsigned int>(natom))
              break;
          }
      }
    else
      {
        avisit.SetBitOn(natom);
        atom = mol.GetAtom(natom);
        for(bond = atom->BeginBond(k);bond;bond=atom->NextBond(k))
          if( !bvisit[bond->GetIdx()])
            {
              path[depth] = bond->GetIdx();
              bvisit.SetBitOn(bond->GetIdx());
              FindRings(mol,path,avisit,bvisit,bond->GetNbrAtomIdx(atom),
                        depth+1);
            }
      }
  }

  void OBRingSearch::RemoveRedundant(int frj)
  {
    OBBitVec tmp;
    int i,j;

    //remove identical rings
    for (i = _rlist.size()-1;i > 0;i--)
      for (j = i-1;j >= 0;j--)
        if ((_rlist[i])->_pathset == (_rlist[j])->_pathset)
          {
            delete _rlist[i];
            _rlist.erase(_rlist.begin()+i);
            break;
          }

    //make sure tmp is the same size as the rings
    for (j = 0;j < (signed)_rlist.size();++j)
      tmp = (_rlist[j])->_pathset;

    //remove larger rings that cover the same atoms as smaller rings
    for (i = _rlist.size()-1;i >= 0;i--)
      {
        tmp.Clear();
        for (j = 0;j < (signed)_rlist.size();++j)
          if ((_rlist[j])->_path.size() <= (_rlist[i])->_path.size() && i != j)
            tmp |= (_rlist[j])->_pathset;

        tmp = tmp & (_rlist[i])->_pathset;

        if (tmp == (_rlist[i])->_pathset)
          {
            delete _rlist[i];
            _rlist.erase(_rlist.begin()+i);
          }

        if (_rlist.size() == (unsigned)frj)
          break;
      }
  }

void testOBRotorListFixedBonds()
{
  // 1 2 3 4 5 6 7 8
  // C-C-C-C-C-C-C-C
  //  0 1 2 3 4 5 6
  OBMolPtr mol = OBTestUtil::ReadFile("octane.cml");

  // test with no bonds fixed
  OBRotorList rlist1;
  rlist1.Setup(*mol);
  OB_ASSERT(rlist1.Size() == 5);

  // test with bond 3 fixed
  OBBitVec fixedBonds;
  fixedBonds.SetBitOn(3);
  rlist1.SetFixedBonds(fixedBonds);
  rlist1.Setup(*mol);
  OB_ASSERT(rlist1.Size() == 4);

  // test with bond 1, 3, 5 fixed
  fixedBonds.SetBitOn(1);
  fixedBonds.SetBitOn(5);
  rlist1.SetFixedBonds(fixedBonds);
  rlist1.Setup(*mol);
  OB_ASSERT(rlist1.Size() == 2);

  // test with bond 1, 2, 3, 5 fixed
  fixedBonds.SetBitOn(2);
  rlist1.SetFixedBonds(fixedBonds);
  rlist1.Setup(*mol);
  OB_ASSERT(rlist1.Size() == 1);




}

void testIsomorphismMask()
{
  // read file: 3 6-rings
  //
  //     /\ /\ /\
  //    |  |  |  |
  //     \/ \/ \/
  //
  OBMol mol;
  OBConversion conv;
  conv.SetInFormat("cml");
  std::ifstream ifs(OBTestUtil::GetFilename("isomorphism1.cml").c_str());
  OB_REQUIRE( ifs );
  conv.Read(&mol, &ifs);

  OBQuery *query = CompileSmilesQuery("C1CCCCC1");
  OBIsomorphismMapper *mapper = OBIsomorphismMapper::GetInstance(query);

  // no mask
  OBIsomorphismMapper::Mappings maps;
  mapper->MapUnique(&mol, maps);
  cout << maps.size() << endl;
  OB_ASSERT( maps.size() == 3 );

  // mask first ring
  OBBitVec mask;
  for (int i = 0; i < 6; ++i)
    mask.SetBitOn(i+1);
  mapper->MapUnique(&mol, maps, mask);
  cout << maps.size() << endl;
  OB_ASSERT( maps.size() == 1 );

  // mask second ring also
  for (int i = 6; i < 10; ++i)
    mask.SetBitOn(i+1);
  mapper->MapUnique(&mol, maps, mask);
  cout << maps.size() << endl;
  OB_ASSERT( maps.size() == 2 );

  // just mask last ring (atomIds 7-8, 10-13)
  mask.Clear();
  for (int i = 10; i < 14; ++i)
    mask.SetBitOn(i+1);
  mask.SetBitOn(7 + 1); mask.SetBitOn(8 + 1);
  mapper->MapUnique(&mol, maps, mask);
  cout << maps.size() << endl;
  OB_ASSERT( maps.size() == 1 ); // Should be same result as masking just the first ring

  delete query;
  delete mapper;
}

  void BuildOBRTreeVector(OBAtom *atom,OBRTree *prv,vector<OBRTree*> &vt,OBBitVec &bv)
  {
    vt[atom->GetIdx()] = new OBRTree (atom,prv);

    int i;
    OBAtom *nbr;
    OBMol *mol = (OBMol*)atom->GetParent();
    OBBitVec curr,used,next;
    vector<OBBond*>::iterator j;
    curr |= atom->GetIdx();
    used = bv|curr;

#define OB_RTREE_CUTOFF 20

    int level=0;
    for (;;)
      {
        next.Clear();
        for (i = curr.NextBit(0);i != bv.EndBit();i = curr.NextBit(i))
          {
            atom = mol->GetAtom(i);
            for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
              if (!used[nbr->GetIdx()])
                {
                  next |= nbr->GetIdx();
                  used |= nbr->GetIdx();
                  vt[nbr->GetIdx()] = new OBRTree (nbr,vt[atom->GetIdx()]);
                }
          }

        if (next.Empty())
          break;
        curr = next;
        level++;
        if (level > OB_RTREE_CUTOFF)
          break;
      }
#undef OB_RTREE_CUTOFF
  }

 /** Return true if \p bv1 and \p bv2 are equivalent
     Not that they may be of different size, and still equivalent provided that the extra bits are all zero.
     \param[in] bv1 A bit vector
     \param[in] bv2 Another bit vector
     \return true if equal, false otherwise
 */
 bool operator== (const OBBitVec & bv1, const OBBitVec & bv2)
 {
   if (bv1.GetSize() < bv2.GetSize())
     { // bv1 smaller than bv2
     unsigned i;
     for (i = 0; i < bv1.GetSize(); ++ i)
       if (bv1._set[i] != bv2._set[i])
         return false;
     for (; i < bv2.GetSize(); ++ i)
       if (bv2._set[i] != 0)
         return false;
     }
   else
     { // bv2 smaller or equal than bv1
     unsigned i;
     for (i = 0; i < bv2.GetSize(); ++ i)
       if (bv1._set[i] != bv2._set[i])
         return false;
     for (; i < bv1.GetSize(); ++ i)
       if (bv1._set[i] != 0)
         return false;
     }
   return true;
 }

  bool OBDepict::DrawMolecule(OBMol *mol)
  {
    if (!d->painter)
      return false;

    d->mol = mol;

    double width=0.0, height=0.0;

    OBAtom *atom;
    OBBondIterator j;
    OBAtomIterator i;

    if(mol->NumAtoms()>0) {
      // scale bond lengths
      double bondLengthSum = 0.0;
      for (OBBond *bond = mol->BeginBond(j); bond; bond = mol->NextBond(j))
        bondLengthSum += bond->GetLength();
      const double averageBondLength = bondLengthSum / mol->NumBonds();
      const double f = mol->NumBonds() ? d->bondLength / averageBondLength : 1.0;
      for (atom = mol->BeginAtom(i); atom; atom = mol->NextAtom(i))
        atom->SetVector(atom->GetX() * f, atom->GetY() * f, 0.0);

      // find min/max values
      double min_x, max_x;
      double min_y, max_y;
      atom = mol->BeginAtom(i);
      min_x = max_x = atom->GetX();
      min_y = max_y = atom->GetY();
      for (atom = mol->NextAtom(i); atom; atom = mol->NextAtom(i)) {
        min_x = std::min(min_x, atom->GetX());
        max_x = std::max(max_x, atom->GetX());
        min_y = std::min(min_y, atom->GetY());
        max_y = std::max(max_y, atom->GetY());
      }

      const double margin = 40.0;
      // translate all atoms so the bottom-left atom is at margin,margin
      for (atom = mol->BeginAtom(i); atom; atom = mol->NextAtom(i))
        atom->SetVector(atom->GetX() - min_x + margin, atom->GetY() - min_y + margin, 0.0);

      width  = max_x - min_x + 2*margin;
      height = max_y - min_y + 2*margin;
      
      //d->painter->SetPenWidth(d->penWidth);
      //d->painter->SetPenColor(d->pen));
      //d->painter->SetFillColor(OBColor("black"));
    }

    d->painter->NewCanvas(width, height);
    
    // draw bonds
    if(d->options & genWedgeHash)
      d->SetWedgeAndHash(mol);
    for (OBBond *bond = mol->BeginBond(j); bond; bond = mol->NextBond(j)) {
      OBAtom *begin = bond->GetBeginAtom();
      OBAtom *end = bond->GetEndAtom();

      if((d->options & internalColor) && bond->HasData("color"))
        d->painter->SetPenColor(OBColor(bond->GetData("color")->GetValue()));
      else
        d->painter->SetPenColor(d->bondColor);

      if (bond->IsWedge()) {
        d->DrawWedge(begin, end);
      } else if (bond->IsHash()) {
        d->DrawHash(begin, end);
      } else if (!bond->IsInRing()) {
        d->DrawSimpleBond(begin, end, bond->GetBO());
      }
    }
    
    // draw ring bonds
    std::vector<OBRing*> rings(mol->GetSSSR());
    OBBitVec drawnBonds;
    for (std::vector<OBRing*>::iterator k = rings.begin(); k != rings.end(); ++k) {
      OBRing *ring = *k;
      std::vector<int> indexes = ring->_path;
      vector3 center(VZero);
      for (std::vector<int>::iterator l = indexes.begin(); l != indexes.end(); ++l) {
        center += mol->GetAtom(*l)->GetVector();        
      }
      center /= indexes.size();

      for (unsigned int l = 0; l < indexes.size(); ++l) {
        OBAtom *begin = mol->GetAtom(indexes[l]);
        OBAtom *end;
        if (l+1 < indexes.size())
          end = mol->GetAtom(indexes[l+1]);
        else
          end = mol->GetAtom(indexes[0]);

        OBBond *ringBond = mol->GetBond(begin, end);
        if (drawnBonds.BitIsSet(ringBond->GetId()))
          continue;

        if((d->options & internalColor) && ringBond->HasData("color"))
          d->painter->SetPenColor(OBColor(ringBond->GetData("color")->GetValue()));
        else
          d->painter->SetPenColor(d->bondColor);
//.........这里部分代码省略.........