[bcf653] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | */
|
---|
| 7 |
|
---|
[cee0b57] | 8 | /*
|
---|
| 9 | * molecule_graph.cpp
|
---|
| 10 | *
|
---|
| 11 | * Created on: Oct 5, 2009
|
---|
| 12 | * Author: heber
|
---|
| 13 | */
|
---|
| 14 |
|
---|
[bf3817] | 15 | // include config.h
|
---|
[aafd77] | 16 | #ifdef HAVE_CONFIG_H
|
---|
| 17 | #include <config.h>
|
---|
| 18 | #endif
|
---|
| 19 |
|
---|
[ad011c] | 20 | #include "CodePatterns/MemDebug.hpp"
|
---|
[112b09] | 21 |
|
---|
[a564be] | 22 | #include <stack>
|
---|
| 23 |
|
---|
[f66195] | 24 | #include "atom.hpp"
|
---|
| 25 | #include "bond.hpp"
|
---|
[b70721] | 26 | #include "bondgraph.hpp"
|
---|
[cee0b57] | 27 | #include "config.hpp"
|
---|
[f66195] | 28 | #include "element.hpp"
|
---|
[1d5afa5] | 29 | #include "Helpers/defs.hpp"
|
---|
| 30 | #include "Helpers/fast_functions.hpp"
|
---|
[952f38] | 31 | #include "Helpers/helpers.hpp"
|
---|
[1d5afa5] | 32 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
[b8b75d] | 33 | #include "linkedcell.hpp"
|
---|
[f66195] | 34 | #include "lists.hpp"
|
---|
[cee0b57] | 35 | #include "molecule.hpp"
|
---|
[b34306] | 36 | #include "World.hpp"
|
---|
[9d83b6] | 37 | #include "WorldTime.hpp"
|
---|
[84c494] | 38 | #include "Box.hpp"
|
---|
[cee0b57] | 39 |
|
---|
[1d5afa5] | 40 | #include "CodePatterns/Assert.hpp"
|
---|
| 41 | #include "CodePatterns/Info.hpp"
|
---|
| 42 | #include "CodePatterns/Log.hpp"
|
---|
| 43 | #include "CodePatterns/Verbose.hpp"
|
---|
| 44 |
|
---|
| 45 | #define MAXBONDS 8
|
---|
| 46 |
|
---|
[9eefda] | 47 | struct BFSAccounting
|
---|
| 48 | {
|
---|
| 49 | atom **PredecessorList;
|
---|
| 50 | int *ShortestPathList;
|
---|
| 51 | enum Shading *ColorList;
|
---|
[a564be] | 52 | std::deque<atom *> *BFSStack;
|
---|
| 53 | std::deque<atom *> *TouchedStack;
|
---|
[9eefda] | 54 | int AtomCount;
|
---|
| 55 | int BondOrder;
|
---|
| 56 | atom *Root;
|
---|
| 57 | bool BackStepping;
|
---|
| 58 | int CurrentGraphNr;
|
---|
| 59 | int ComponentNr;
|
---|
| 60 | };
|
---|
[cee0b57] | 61 |
|
---|
[9eefda] | 62 | /** Accounting data for Depth First Search.
|
---|
| 63 | */
|
---|
| 64 | struct DFSAccounting
|
---|
| 65 | {
|
---|
[a564be] | 66 | std::deque<atom *> *AtomStack;
|
---|
| 67 | std::deque<bond *> *BackEdgeStack;
|
---|
[9eefda] | 68 | int CurrentGraphNr;
|
---|
| 69 | int ComponentNumber;
|
---|
| 70 | atom *Root;
|
---|
| 71 | bool BackStepping;
|
---|
| 72 | };
|
---|
| 73 |
|
---|
| 74 | /************************************* Functions for class molecule *********************************/
|
---|
[cee0b57] | 75 |
|
---|
| 76 | /** Creates an adjacency list of the molecule.
|
---|
| 77 | * We obtain an outside file with the indices of atoms which are bondmembers.
|
---|
| 78 | */
|
---|
[e138de] | 79 | void molecule::CreateAdjacencyListFromDbondFile(ifstream *input)
|
---|
[cee0b57] | 80 | {
|
---|
[c68c90] | 81 | Info FunctionInfo(__func__);
|
---|
[cee0b57] | 82 | // 1 We will parse bonds out of the dbond file created by tremolo.
|
---|
[44a59b] | 83 | int atom1, atom2;
|
---|
| 84 | atom *Walker, *OtherWalker;
|
---|
[c68c90] | 85 | char line[MAXSTRINGSIZE];
|
---|
[44a59b] | 86 |
|
---|
[c68c90] | 87 | if (input->fail()) {
|
---|
| 88 | DoeLog(0) && (eLog() << Verbose(0) << "Opening of bond file failed \n");
|
---|
| 89 | performCriticalExit();
|
---|
[44a59b] | 90 | };
|
---|
[bd6bfa] | 91 | doCountAtoms();
|
---|
[44a59b] | 92 |
|
---|
[c68c90] | 93 | // skip header
|
---|
| 94 | input->getline(line,MAXSTRINGSIZE);
|
---|
| 95 | DoLog(1) && (Log() << Verbose(1) << "Scanning file ... \n");
|
---|
[44a59b] | 96 | while (!input->eof()) // Check whether we read everything already
|
---|
| 97 | {
|
---|
[c68c90] | 98 | input->getline(line,MAXSTRINGSIZE);
|
---|
| 99 | stringstream zeile(line);
|
---|
| 100 | zeile >> atom1;
|
---|
| 101 | zeile >> atom2;
|
---|
[44a59b] | 102 |
|
---|
[c68c90] | 103 | DoLog(2) && (Log() << Verbose(2) << "Looking for atoms " << atom1 << " and " << atom2 << "." << endl);
|
---|
[9eefda] | 104 | if (atom2 < atom1) //Sort indices of atoms in order
|
---|
[a0064e] | 105 | std::swap(atom1, atom2);
|
---|
[9eefda] | 106 | Walker = FindAtom(atom1);
|
---|
[05a97c] | 107 | ASSERT(Walker,"Could not find an atom with the ID given in dbond file");
|
---|
[9eefda] | 108 | OtherWalker = FindAtom(atom2);
|
---|
[05a97c] | 109 | ASSERT(OtherWalker,"Could not find an atom with the ID given in dbond file");
|
---|
[44a59b] | 110 | AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
|
---|
| 111 | }
|
---|
[9eefda] | 112 | }
|
---|
| 113 | ;
|
---|
[cee0b57] | 114 |
|
---|
| 115 | /** Creates an adjacency list of the molecule.
|
---|
| 116 | * Generally, we use the CSD approach to bond recognition, that is the the distance
|
---|
| 117 | * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
|
---|
| 118 | * a threshold t = 0.4 Angstroem.
|
---|
| 119 | * To make it O(N log N) the function uses the linked-cell technique as follows:
|
---|
| 120 | * The procedure is step-wise:
|
---|
| 121 | * -# Remove every bond in list
|
---|
| 122 | * -# Count the atoms in the molecule with CountAtoms()
|
---|
| 123 | * -# partition cell into smaller linked cells of size \a bonddistance
|
---|
| 124 | * -# put each atom into its corresponding cell
|
---|
| 125 | * -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
|
---|
| 126 | * -# correct the bond degree iteratively (single->double->triple bond)
|
---|
| 127 | * -# finally print the bond list to \a *out if desired
|
---|
| 128 | * \param *out out stream for printing the matrix, NULL if no output
|
---|
| 129 | * \param bonddistance length of linked cells (i.e. maximum minimal length checked)
|
---|
| 130 | * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
|
---|
[b70721] | 131 | * \param *minmaxdistance function to give upper and lower bound on whether particle is bonded to some other
|
---|
| 132 | * \param *BG BondGraph with the member function above or NULL, if just standard covalent should be used.
|
---|
[cee0b57] | 133 | */
|
---|
[e138de] | 134 | void molecule::CreateAdjacencyList(double bonddistance, bool IsAngstroem, void (BondGraph::*minmaxdistance)(BondedParticle * const , BondedParticle * const , double &, double &, bool), BondGraph *BG)
|
---|
[cee0b57] | 135 | {
|
---|
[b8b75d] | 136 | atom *Walker = NULL;
|
---|
| 137 | atom *OtherWalker = NULL;
|
---|
| 138 | int n[NDIM];
|
---|
[b70721] | 139 | double MinDistance, MaxDistance;
|
---|
[b8b75d] | 140 | LinkedCell *LC = NULL;
|
---|
[b70721] | 141 | bool free_BG = false;
|
---|
[014475] | 142 | Box &domain = World::getInstance().getDomain();
|
---|
[b70721] | 143 |
|
---|
| 144 | if (BG == NULL) {
|
---|
| 145 | BG = new BondGraph(IsAngstroem);
|
---|
| 146 | free_BG = true;
|
---|
| 147 | }
|
---|
[cee0b57] | 148 |
|
---|
| 149 | BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
|
---|
[a67d19] | 150 | DoLog(0) && (Log() << Verbose(0) << "Begin of CreateAdjacencyList." << endl);
|
---|
[cee0b57] | 151 | // remove every bond from the list
|
---|
[9d83b6] | 152 | for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 153 | BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 154 | for(BondList::iterator BondRunner = ListOfBonds.begin();
|
---|
| 155 | !ListOfBonds.empty();
|
---|
| 156 | BondRunner = ListOfBonds.begin())
|
---|
[e08c46] | 157 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
| 158 | delete((*BondRunner));
|
---|
[9d83b6] | 159 | }
|
---|
[3c349b] | 160 | BondCount = 0;
|
---|
[cee0b57] | 161 |
|
---|
| 162 | // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
|
---|
[a7b761b] | 163 | DoLog(1) && (Log() << Verbose(1) << "AtomCount " << getAtomCount() << " and bonddistance is " << bonddistance << "." << endl);
|
---|
[cee0b57] | 164 |
|
---|
[c66537] | 165 | if ((getAtomCount() > 1) && (bonddistance > 0.1)) {
|
---|
[a67d19] | 166 | DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
|
---|
[af2c424] | 167 | LC = new LinkedCell(*this, bonddistance);
|
---|
[cee0b57] | 168 |
|
---|
[b8b75d] | 169 | // create a list to map Tesselpoint::nr to atom *
|
---|
[a67d19] | 170 | DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl);
|
---|
[f2bb0f] | 171 |
|
---|
[53731f] | 172 | // set numbers for atoms that can later be used
|
---|
| 173 | int i=0;
|
---|
| 174 | for(internal_iterator iter = atoms.begin();iter!= atoms.end(); ++iter){
|
---|
| 175 | (*iter)->nr = i++;
|
---|
[cee0b57] | 176 | }
|
---|
| 177 |
|
---|
| 178 | // 3a. go through every cell
|
---|
[a67d19] | 179 | DoLog(2) && (Log() << Verbose(2) << "Celling ... " << endl);
|
---|
[b8b75d] | 180 | for (LC->n[0] = 0; LC->n[0] < LC->N[0]; LC->n[0]++)
|
---|
| 181 | for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
|
---|
| 182 | for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
|
---|
[734816] | 183 | const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
|
---|
[4e855e] | 184 | Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
|
---|
[b8b75d] | 185 | if (List != NULL) {
|
---|
[734816] | 186 | for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
---|
[f2bb0f] | 187 | Walker = dynamic_cast<atom*>(*Runner);
|
---|
| 188 | ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
|
---|
[4e855e] | 189 | Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
|
---|
[cee0b57] | 190 | // 3c. check for possible bond between each atom in this and every one in the 27 cells
|
---|
[9eefda] | 191 | for (n[0] = -1; n[0] <= 1; n[0]++)
|
---|
| 192 | for (n[1] = -1; n[1] <= 1; n[1]++)
|
---|
| 193 | for (n[2] = -1; n[2] <= 1; n[2]++) {
|
---|
[734816] | 194 | const LinkedCell::LinkedNodes *OtherList = LC->GetRelativeToCurrentCell(n);
|
---|
[b8b75d] | 195 | if (OtherList != NULL) {
|
---|
[4e855e] | 196 | Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
|
---|
[734816] | 197 | for (LinkedCell::LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
|
---|
[b8b75d] | 198 | if ((*OtherRunner)->nr > Walker->nr) {
|
---|
[f2bb0f] | 199 | OtherWalker = dynamic_cast<atom*>(*OtherRunner);
|
---|
| 200 | ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
|
---|
[e5ad5c] | 201 | (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
|
---|
[d74077] | 202 | const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
|
---|
[4e855e] | 203 | Log() << Verbose(1) << "Checking distance " << distance << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
|
---|
[b70721] | 204 | const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
|
---|
[4e855e] | 205 | Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
|
---|
[e5ad5c] | 206 | if (OtherWalker->father->nr > Walker->father->nr) {
|
---|
| 207 | if (status) { // create bond if distance is smaller
|
---|
[4e855e] | 208 | Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
|
---|
[e5ad5c] | 209 | AddBond(Walker->father, OtherWalker->father, 1); // also increases molecule::BondCount
|
---|
| 210 | } else {
|
---|
[4e855e] | 211 | Log() << Verbose(1) << "Not Adding: distance too great." << endl;
|
---|
[e5ad5c] | 212 | }
|
---|
[b8b75d] | 213 | } else {
|
---|
[4e855e] | 214 | Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
|
---|
[b8b75d] | 215 | }
|
---|
[cee0b57] | 216 | }
|
---|
| 217 | }
|
---|
| 218 | }
|
---|
| 219 | }
|
---|
| 220 | }
|
---|
| 221 | }
|
---|
| 222 | }
|
---|
[9eefda] | 223 | delete (LC);
|
---|
[a67d19] | 224 | DoLog(1) && (Log() << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl);
|
---|
[cee0b57] | 225 |
|
---|
[b8b75d] | 226 | // correct bond degree by comparing valence and bond degree
|
---|
[a67d19] | 227 | DoLog(2) && (Log() << Verbose(2) << "Correcting bond degree ... " << endl);
|
---|
[e138de] | 228 | CorrectBondDegree();
|
---|
[cee0b57] | 229 |
|
---|
[b8b75d] | 230 | // output bonds for debugging (if bond chain list was correctly installed)
|
---|
[c743f8] | 231 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
|
---|
[b8b75d] | 232 | } else
|
---|
[a7b761b] | 233 | DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << getAtomCount() << ", thus no bonds, no connections!." << endl);
|
---|
[a67d19] | 234 | DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl);
|
---|
[b70721] | 235 | if (free_BG)
|
---|
| 236 | delete(BG);
|
---|
[9eefda] | 237 | }
|
---|
| 238 | ;
|
---|
[cee0b57] | 239 |
|
---|
[e08c46] | 240 | /** Checks for presence of bonds within atom list.
|
---|
| 241 | * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
|
---|
| 242 | * \return true - bonds present, false - no bonds
|
---|
| 243 | */
|
---|
[e4afb4] | 244 | bool molecule::hasBondStructure() const
|
---|
[e08c46] | 245 | {
|
---|
[9d83b6] | 246 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 247 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 248 | if (!ListOfBonds.empty())
|
---|
[e08c46] | 249 | return true;
|
---|
[9d83b6] | 250 | }
|
---|
[e08c46] | 251 | return false;
|
---|
| 252 | }
|
---|
| 253 |
|
---|
| 254 | /** Counts the number of present bonds.
|
---|
| 255 | * \return number of bonds
|
---|
| 256 | */
|
---|
| 257 | unsigned int molecule::CountBonds() const
|
---|
| 258 | {
|
---|
| 259 | unsigned int counter = 0;
|
---|
[9d83b6] | 260 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 261 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 262 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 263 | BondRunner != ListOfBonds.end();
|
---|
| 264 | ++BondRunner)
|
---|
[e08c46] | 265 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
| 266 | counter++;
|
---|
[9d83b6] | 267 | }
|
---|
[e08c46] | 268 | return counter;
|
---|
| 269 | }
|
---|
| 270 |
|
---|
[b8b75d] | 271 | /** Prints a list of all bonds to \a *out.
|
---|
| 272 | * \param output stream
|
---|
| 273 | */
|
---|
[e138de] | 274 | void molecule::OutputBondsList() const
|
---|
[b8b75d] | 275 | {
|
---|
[a67d19] | 276 | DoLog(1) && (Log() << Verbose(1) << endl << "From contents of bond chain list:");
|
---|
[9d83b6] | 277 | for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) {
|
---|
| 278 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 279 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 280 | BondRunner != ListOfBonds.end();
|
---|
| 281 | ++BondRunner)
|
---|
[e08c46] | 282 | if ((*BondRunner)->leftatom == *AtomRunner) {
|
---|
| 283 | DoLog(0) && (Log() << Verbose(0) << *(*BondRunner) << "\t" << endl);
|
---|
| 284 | }
|
---|
[9d83b6] | 285 | }
|
---|
[a67d19] | 286 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
[9eefda] | 287 | }
|
---|
| 288 | ;
|
---|
[cee0b57] | 289 |
|
---|
[b8b75d] | 290 | /** correct bond degree by comparing valence and bond degree.
|
---|
| 291 | * correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
|
---|
| 292 | * iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
|
---|
| 293 | * preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
|
---|
| 294 | * double bonds as was expected.
|
---|
| 295 | * \param *out output stream for debugging
|
---|
| 296 | * \return number of bonds that could not be corrected
|
---|
| 297 | */
|
---|
[e138de] | 298 | int molecule::CorrectBondDegree() const
|
---|
[b8b75d] | 299 | {
|
---|
[99593f] | 300 | int No = 0, OldNo = -1;
|
---|
[b8b75d] | 301 |
|
---|
| 302 | if (BondCount != 0) {
|
---|
[a67d19] | 303 | DoLog(1) && (Log() << Verbose(1) << "Correcting Bond degree of each bond ... " << endl);
|
---|
[b8b75d] | 304 | do {
|
---|
[99593f] | 305 | OldNo = No;
|
---|
[00ef5c] | 306 | No=0;
|
---|
| 307 | BOOST_FOREACH(atom *atom,atoms){
|
---|
| 308 | No+=atom->CorrectBondDegree();
|
---|
| 309 | }
|
---|
[99593f] | 310 | } while (OldNo != No);
|
---|
[a67d19] | 311 | DoLog(0) && (Log() << Verbose(0) << " done." << endl);
|
---|
[b8b75d] | 312 | } else {
|
---|
[a7b761b] | 313 | DoLog(1) && (Log() << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << getAtomCount() << " atoms." << endl);
|
---|
[b8b75d] | 314 | }
|
---|
[a67d19] | 315 | DoLog(0) && (Log() << Verbose(0) << No << " bonds could not be corrected." << endl);
|
---|
[cee0b57] | 316 |
|
---|
[266237] | 317 | return (No);
|
---|
[9eefda] | 318 | }
|
---|
| 319 | ;
|
---|
[cee0b57] | 320 |
|
---|
| 321 | /** Counts all cyclic bonds and returns their number.
|
---|
| 322 | * \note Hydrogen bonds can never by cyclic, thus no check for that
|
---|
| 323 | * \param *out output stream for debugging
|
---|
| 324 | * \return number opf cyclic bonds
|
---|
| 325 | */
|
---|
[e138de] | 326 | int molecule::CountCyclicBonds()
|
---|
[cee0b57] | 327 | {
|
---|
[266237] | 328 | NoCyclicBonds = 0;
|
---|
[cee0b57] | 329 | int *MinimumRingSize = NULL;
|
---|
| 330 | MoleculeLeafClass *Subgraphs = NULL;
|
---|
[a564be] | 331 | std::deque<bond *> *BackEdgeStack = NULL;
|
---|
[9d83b6] | 332 | for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 333 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 334 | if ((!ListOfBonds.empty()) && ((*ListOfBonds.begin())->Type == Undetermined)) {
|
---|
[e08c46] | 335 | DoLog(0) && (Log() << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl);
|
---|
| 336 | Subgraphs = DepthFirstSearchAnalysis(BackEdgeStack);
|
---|
| 337 | while (Subgraphs->next != NULL) {
|
---|
| 338 | Subgraphs = Subgraphs->next;
|
---|
| 339 | delete (Subgraphs->previous);
|
---|
| 340 | }
|
---|
| 341 | delete (Subgraphs);
|
---|
| 342 | delete[] (MinimumRingSize);
|
---|
| 343 | break;
|
---|
[cee0b57] | 344 | }
|
---|
[9d83b6] | 345 | }
|
---|
| 346 | for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 347 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 348 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 349 | BondRunner != ListOfBonds.end();
|
---|
| 350 | ++BondRunner)
|
---|
[e08c46] | 351 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
| 352 | if ((*BondRunner)->Cyclic)
|
---|
| 353 | NoCyclicBonds++;
|
---|
[9d83b6] | 354 | }
|
---|
[9eefda] | 355 | delete (BackEdgeStack);
|
---|
[266237] | 356 | return NoCyclicBonds;
|
---|
[9eefda] | 357 | }
|
---|
| 358 | ;
|
---|
[b8b75d] | 359 |
|
---|
[cee0b57] | 360 | /** Returns Shading as a char string.
|
---|
| 361 | * \param color the Shading
|
---|
| 362 | * \return string of the flag
|
---|
| 363 | */
|
---|
[fa649a] | 364 | string molecule::GetColor(enum Shading color) const
|
---|
[cee0b57] | 365 | {
|
---|
[9eefda] | 366 | switch (color) {
|
---|
[cee0b57] | 367 | case white:
|
---|
| 368 | return "white";
|
---|
| 369 | break;
|
---|
| 370 | case lightgray:
|
---|
| 371 | return "lightgray";
|
---|
| 372 | break;
|
---|
| 373 | case darkgray:
|
---|
| 374 | return "darkgray";
|
---|
| 375 | break;
|
---|
| 376 | case black:
|
---|
| 377 | return "black";
|
---|
| 378 | break;
|
---|
| 379 | default:
|
---|
| 380 | return "uncolored";
|
---|
| 381 | break;
|
---|
| 382 | };
|
---|
[9eefda] | 383 | }
|
---|
| 384 | ;
|
---|
[cee0b57] | 385 |
|
---|
[9eefda] | 386 | /** Sets atom::GraphNr and atom::LowpointNr to BFSAccounting::CurrentGraphNr.
|
---|
| 387 | * \param *out output stream for debugging
|
---|
| 388 | * \param *Walker current node
|
---|
| 389 | * \param &BFS structure with accounting data for BFS
|
---|
| 390 | */
|
---|
[e138de] | 391 | void DepthFirstSearchAnalysis_SetWalkersGraphNr(atom *&Walker, struct DFSAccounting &DFS)
|
---|
[174e0e] | 392 | {
|
---|
[9eefda] | 393 | if (!DFS.BackStepping) { // if we don't just return from (8)
|
---|
| 394 | Walker->GraphNr = DFS.CurrentGraphNr;
|
---|
| 395 | Walker->LowpointNr = DFS.CurrentGraphNr;
|
---|
[68f03d] | 396 | DoLog(1) && (Log() << Verbose(1) << "Setting Walker[" << Walker->getName() << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "." << endl);
|
---|
[a564be] | 397 | DFS.AtomStack->push_front(Walker);
|
---|
[9eefda] | 398 | DFS.CurrentGraphNr++;
|
---|
[174e0e] | 399 | }
|
---|
[9eefda] | 400 | }
|
---|
| 401 | ;
|
---|
[174e0e] | 402 |
|
---|
[9eefda] | 403 | /** During DFS goes along unvisited bond and touches other atom.
|
---|
| 404 | * Sets bond::type, if
|
---|
| 405 | * -# BackEdge: set atom::LowpointNr and push on \a BackEdgeStack
|
---|
| 406 | * -# TreeEgde: set atom::Ancestor and continue with Walker along this edge
|
---|
| 407 | * Continue until molecule::FindNextUnused() finds no more unused bonds.
|
---|
| 408 | * \param *out output stream for debugging
|
---|
| 409 | * \param *mol molecule with atoms and finding unused bonds
|
---|
| 410 | * \param *&Binder current edge
|
---|
| 411 | * \param &DFS DFS accounting data
|
---|
| 412 | */
|
---|
[e138de] | 413 | void DepthFirstSearchAnalysis_ProbeAlongUnusedBond(const molecule * const mol, atom *&Walker, bond *&Binder, struct DFSAccounting &DFS)
|
---|
[174e0e] | 414 | {
|
---|
| 415 | atom *OtherAtom = NULL;
|
---|
| 416 |
|
---|
| 417 | do { // (3) if Walker has no unused egdes, go to (5)
|
---|
[9eefda] | 418 | DFS.BackStepping = false; // reset backstepping flag for (8)
|
---|
[174e0e] | 419 | if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
|
---|
| 420 | Binder = mol->FindNextUnused(Walker);
|
---|
| 421 | if (Binder == NULL)
|
---|
| 422 | break;
|
---|
[a67d19] | 423 | DoLog(2) && (Log() << Verbose(2) << "Current Unused Bond is " << *Binder << "." << endl);
|
---|
[174e0e] | 424 | // (4) Mark Binder used, ...
|
---|
| 425 | Binder->MarkUsed(black);
|
---|
| 426 | OtherAtom = Binder->GetOtherAtom(Walker);
|
---|
[68f03d] | 427 | DoLog(2) && (Log() << Verbose(2) << "(4) OtherAtom is " << OtherAtom->getName() << "." << endl);
|
---|
[174e0e] | 428 | if (OtherAtom->GraphNr != -1) {
|
---|
| 429 | // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
|
---|
| 430 | Binder->Type = BackEdge;
|
---|
[a564be] | 431 | DFS.BackEdgeStack->push_front(Binder);
|
---|
[9eefda] | 432 | Walker->LowpointNr = (Walker->LowpointNr < OtherAtom->GraphNr) ? Walker->LowpointNr : OtherAtom->GraphNr;
|
---|
[68f03d] | 433 | DoLog(3) && (Log() << Verbose(3) << "(4a) Visited: Setting Lowpoint of Walker[" << Walker->getName() << "] to " << Walker->LowpointNr << "." << endl);
|
---|
[174e0e] | 434 | } else {
|
---|
| 435 | // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
|
---|
| 436 | Binder->Type = TreeEdge;
|
---|
| 437 | OtherAtom->Ancestor = Walker;
|
---|
| 438 | Walker = OtherAtom;
|
---|
[68f03d] | 439 | DoLog(3) && (Log() << Verbose(3) << "(4b) Not Visited: OtherAtom[" << OtherAtom->getName() << "]'s Ancestor is now " << OtherAtom->Ancestor->getName() << ", Walker is OtherAtom " << OtherAtom->getName() << "." << endl);
|
---|
[174e0e] | 440 | break;
|
---|
| 441 | }
|
---|
| 442 | Binder = NULL;
|
---|
[9eefda] | 443 | } while (1); // (3)
|
---|
| 444 | }
|
---|
| 445 | ;
|
---|
[174e0e] | 446 |
|
---|
[9eefda] | 447 | /** Checks whether we have a new component.
|
---|
| 448 | * if atom::LowpointNr of \a *&Walker is greater than atom::GraphNr of its atom::Ancestor, we have a new component.
|
---|
| 449 | * Meaning that if we touch upon a node who suddenly has a smaller atom::LowpointNr than its ancestor, then we
|
---|
| 450 | * have a found a new branch in the graph tree.
|
---|
| 451 | * \param *out output stream for debugging
|
---|
| 452 | * \param *mol molecule with atoms and finding unused bonds
|
---|
| 453 | * \param *&Walker current node
|
---|
| 454 | * \param &DFS DFS accounting data
|
---|
| 455 | */
|
---|
[e138de] | 456 | void DepthFirstSearchAnalysis_CheckForaNewComponent(const molecule * const mol, atom *&Walker, struct DFSAccounting &DFS, MoleculeLeafClass *&LeafWalker)
|
---|
[174e0e] | 457 | {
|
---|
| 458 | atom *OtherAtom = NULL;
|
---|
| 459 |
|
---|
| 460 | // (5) if Ancestor of Walker is ...
|
---|
[68f03d] | 461 | DoLog(1) && (Log() << Verbose(1) << "(5) Number of Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "] is " << Walker->Ancestor->GraphNr << "." << endl);
|
---|
[174e0e] | 462 |
|
---|
[9eefda] | 463 | if (Walker->Ancestor->GraphNr != DFS.Root->GraphNr) {
|
---|
[174e0e] | 464 | // (6) (Ancestor of Walker is not Root)
|
---|
| 465 | if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
|
---|
| 466 | // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
|
---|
| 467 | Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
|
---|
[68f03d] | 468 | DoLog(2) && (Log() << Verbose(2) << "(6) Setting Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "." << endl);
|
---|
[174e0e] | 469 | } else {
|
---|
| 470 | // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
|
---|
| 471 | Walker->Ancestor->SeparationVertex = true;
|
---|
[68f03d] | 472 | DoLog(2) && (Log() << Verbose(2) << "(7) Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s is a separating vertex, creating component." << endl);
|
---|
[9eefda] | 473 | mol->SetNextComponentNumber(Walker->Ancestor, DFS.ComponentNumber);
|
---|
[68f03d] | 474 | DoLog(3) && (Log() << Verbose(3) << "(7) Walker[" << Walker->getName() << "]'s Ancestor's Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
[9eefda] | 475 | mol->SetNextComponentNumber(Walker, DFS.ComponentNumber);
|
---|
[68f03d] | 476 | DoLog(3) && (Log() << Verbose(3) << "(7) Walker[" << Walker->getName() << "]'s Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
[174e0e] | 477 | do {
|
---|
[a564be] | 478 | ASSERT(!DFS.AtomStack->empty(), "DepthFirstSearchAnalysis_CheckForaNewComponent() - DFS.AtomStack is empty!");
|
---|
| 479 | OtherAtom = DFS.AtomStack->front();
|
---|
| 480 | DFS.AtomStack->pop_front();
|
---|
[174e0e] | 481 | LeafWalker->Leaf->AddCopyAtom(OtherAtom);
|
---|
[9eefda] | 482 | mol->SetNextComponentNumber(OtherAtom, DFS.ComponentNumber);
|
---|
[68f03d] | 483 | DoLog(3) && (Log() << Verbose(3) << "(7) Other[" << OtherAtom->getName() << "]'s Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
[174e0e] | 484 | } while (OtherAtom != Walker);
|
---|
[9eefda] | 485 | DFS.ComponentNumber++;
|
---|
[174e0e] | 486 | }
|
---|
| 487 | // (8) Walker becomes its Ancestor, go to (3)
|
---|
[68f03d] | 488 | DoLog(2) && (Log() << Verbose(2) << "(8) Walker[" << Walker->getName() << "] is now its Ancestor " << Walker->Ancestor->getName() << ", backstepping. " << endl);
|
---|
[174e0e] | 489 | Walker = Walker->Ancestor;
|
---|
[9eefda] | 490 | DFS.BackStepping = true;
|
---|
[174e0e] | 491 | }
|
---|
[9eefda] | 492 | }
|
---|
| 493 | ;
|
---|
[174e0e] | 494 |
|
---|
[9eefda] | 495 | /** Cleans the root stack when we have found a component.
|
---|
| 496 | * If we are not DFSAccounting::BackStepping, then we clear the root stack by putting everything into a
|
---|
| 497 | * component down till we meet DFSAccounting::Root.
|
---|
| 498 | * \param *out output stream for debugging
|
---|
| 499 | * \param *mol molecule with atoms and finding unused bonds
|
---|
| 500 | * \param *&Walker current node
|
---|
| 501 | * \param *&Binder current edge
|
---|
| 502 | * \param &DFS DFS accounting data
|
---|
| 503 | */
|
---|
[e138de] | 504 | void DepthFirstSearchAnalysis_CleanRootStackDownTillWalker(const molecule * const mol, atom *&Walker, bond *&Binder, struct DFSAccounting &DFS, MoleculeLeafClass *&LeafWalker)
|
---|
[174e0e] | 505 | {
|
---|
| 506 | atom *OtherAtom = NULL;
|
---|
| 507 |
|
---|
[9eefda] | 508 | if (!DFS.BackStepping) { // coming from (8) want to go to (3)
|
---|
[174e0e] | 509 | // (9) remove all from stack till Walker (including), these and Root form a component
|
---|
[99593f] | 510 | //DFS.AtomStack->Output(out);
|
---|
[9eefda] | 511 | mol->SetNextComponentNumber(DFS.Root, DFS.ComponentNumber);
|
---|
[68f03d] | 512 | DoLog(3) && (Log() << Verbose(3) << "(9) Root[" << DFS.Root->getName() << "]'s Component is " << DFS.ComponentNumber << "." << endl);
|
---|
[9eefda] | 513 | mol->SetNextComponentNumber(Walker, DFS.ComponentNumber);
|
---|
[68f03d] | 514 | DoLog(3) && (Log() << Verbose(3) << "(9) Walker[" << Walker->getName() << "]'s Component is " << DFS.ComponentNumber << "." << endl);
|
---|
[174e0e] | 515 | do {
|
---|
[a564be] | 516 | ASSERT(!DFS.AtomStack->empty(), "DepthFirstSearchAnalysis_CleanRootStackDownTillWalker() - DFS.AtomStack is empty!");
|
---|
| 517 | OtherAtom = DFS.AtomStack->front();
|
---|
| 518 | DFS.AtomStack->pop_front();
|
---|
[174e0e] | 519 | LeafWalker->Leaf->AddCopyAtom(OtherAtom);
|
---|
[9eefda] | 520 | mol->SetNextComponentNumber(OtherAtom, DFS.ComponentNumber);
|
---|
[a564be] | 521 | DoLog(3) && (Log() << Verbose(3) << "(7) Other[" << OtherAtom->getName() << "]'s Component is " << DFS.ComponentNumber << "." << endl);
|
---|
[174e0e] | 522 | } while (OtherAtom != Walker);
|
---|
[9eefda] | 523 | DFS.ComponentNumber++;
|
---|
[174e0e] | 524 |
|
---|
| 525 | // (11) Root is separation vertex, set Walker to Root and go to (4)
|
---|
[9eefda] | 526 | Walker = DFS.Root;
|
---|
[174e0e] | 527 | Binder = mol->FindNextUnused(Walker);
|
---|
[68f03d] | 528 | DoLog(1) && (Log() << Verbose(1) << "(10) Walker is Root[" << DFS.Root->getName() << "], next Unused Bond is " << Binder << "." << endl);
|
---|
[174e0e] | 529 | if (Binder != NULL) { // Root is separation vertex
|
---|
[a67d19] | 530 | DoLog(1) && (Log() << Verbose(1) << "(11) Root is a separation vertex." << endl);
|
---|
[174e0e] | 531 | Walker->SeparationVertex = true;
|
---|
| 532 | }
|
---|
| 533 | }
|
---|
[9eefda] | 534 | }
|
---|
| 535 | ;
|
---|
| 536 |
|
---|
| 537 | /** Initializes DFSAccounting structure.
|
---|
| 538 | * \param *out output stream for debugging
|
---|
| 539 | * \param &DFS accounting structure to allocate
|
---|
[7218f8] | 540 | * \param *mol molecule with AtomCount, BondCount and all atoms
|
---|
[9eefda] | 541 | */
|
---|
[e138de] | 542 | void DepthFirstSearchAnalysis_Init(struct DFSAccounting &DFS, const molecule * const mol)
|
---|
[9eefda] | 543 | {
|
---|
[a564be] | 544 | DFS.AtomStack = new std::deque<atom *> (mol->getAtomCount());
|
---|
[9eefda] | 545 | DFS.CurrentGraphNr = 0;
|
---|
| 546 | DFS.ComponentNumber = 0;
|
---|
| 547 | DFS.BackStepping = false;
|
---|
[7218f8] | 548 | mol->ResetAllBondsToUnused();
|
---|
[a564be] | 549 | DFS.BackEdgeStack->clear();
|
---|
[9eefda] | 550 | }
|
---|
| 551 | ;
|
---|
[174e0e] | 552 |
|
---|
[9eefda] | 553 | /** Free's DFSAccounting structure.
|
---|
| 554 | * \param *out output stream for debugging
|
---|
| 555 | * \param &DFS accounting structure to free
|
---|
| 556 | */
|
---|
[e138de] | 557 | void DepthFirstSearchAnalysis_Finalize(struct DFSAccounting &DFS)
|
---|
[9eefda] | 558 | {
|
---|
| 559 | delete (DFS.AtomStack);
|
---|
[7218f8] | 560 | // delete (DFS.BackEdgeStack); // DON'T free, see DepthFirstSearchAnalysis(), is returned as allocated
|
---|
[9eefda] | 561 | }
|
---|
| 562 | ;
|
---|
[174e0e] | 563 |
|
---|
[00ef5c] | 564 | void molecule::init_DFS(struct DFSAccounting &DFS) const{
|
---|
| 565 | DepthFirstSearchAnalysis_Init(DFS, this);
|
---|
| 566 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::resetGraphNr));
|
---|
| 567 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::InitComponentNr));
|
---|
| 568 | }
|
---|
| 569 |
|
---|
[cee0b57] | 570 | /** Performs a Depth-First search on this molecule.
|
---|
| 571 | * Marks bonds in molecule as cyclic, bridge, ... and atoms as
|
---|
| 572 | * articulations points, ...
|
---|
| 573 | * We use the algorithm from [Even, Graph Algorithms, p.62].
|
---|
| 574 | * \param *out output stream for debugging
|
---|
[a564be] | 575 | * \param *&BackEdgeStack NULL pointer to std::deque<bond *> with all the found back edges, allocated and filled on return
|
---|
[cee0b57] | 576 | * \return list of each disconnected subgraph as an individual molecule class structure
|
---|
| 577 | */
|
---|
[a564be] | 578 | MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(std::deque<bond *> *&BackEdgeStack) const
|
---|
[cee0b57] | 579 | {
|
---|
[9eefda] | 580 | struct DFSAccounting DFS;
|
---|
[a564be] | 581 | BackEdgeStack = new std::deque<bond *> (BondCount);
|
---|
[9eefda] | 582 | DFS.BackEdgeStack = BackEdgeStack;
|
---|
[cee0b57] | 583 | MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
|
---|
| 584 | MoleculeLeafClass *LeafWalker = SubGraphs;
|
---|
[9eefda] | 585 | int OldGraphNr = 0;
|
---|
[174e0e] | 586 | atom *Walker = NULL;
|
---|
[cee0b57] | 587 | bond *Binder = NULL;
|
---|
| 588 |
|
---|
[a7b761b] | 589 | if (getAtomCount() == 0)
|
---|
[046783] | 590 | return SubGraphs;
|
---|
[a67d19] | 591 | DoLog(0) && (Log() << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl);
|
---|
[00ef5c] | 592 | init_DFS(DFS);
|
---|
[cee0b57] | 593 |
|
---|
[9879f6] | 594 | for (molecule::const_iterator iter = begin(); iter != end();) {
|
---|
| 595 | DFS.Root = *iter;
|
---|
[7218f8] | 596 | // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all
|
---|
[a564be] | 597 | DFS.AtomStack->clear();
|
---|
[cee0b57] | 598 |
|
---|
| 599 | // put into new subgraph molecule and add this to list of subgraphs
|
---|
| 600 | LeafWalker = new MoleculeLeafClass(LeafWalker);
|
---|
[5f612ee] | 601 | LeafWalker->Leaf = World::getInstance().createMolecule();
|
---|
[9eefda] | 602 | LeafWalker->Leaf->AddCopyAtom(DFS.Root);
|
---|
[cee0b57] | 603 |
|
---|
[9eefda] | 604 | OldGraphNr = DFS.CurrentGraphNr;
|
---|
| 605 | Walker = DFS.Root;
|
---|
[cee0b57] | 606 | do { // (10)
|
---|
| 607 | do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom
|
---|
[e138de] | 608 | DepthFirstSearchAnalysis_SetWalkersGraphNr(Walker, DFS);
|
---|
[174e0e] | 609 |
|
---|
[e138de] | 610 | DepthFirstSearchAnalysis_ProbeAlongUnusedBond(this, Walker, Binder, DFS);
|
---|
[174e0e] | 611 |
|
---|
[cee0b57] | 612 | if (Binder == NULL) {
|
---|
[a67d19] | 613 | DoLog(2) && (Log() << Verbose(2) << "No more Unused Bonds." << endl);
|
---|
[cee0b57] | 614 | break;
|
---|
| 615 | } else
|
---|
| 616 | Binder = NULL;
|
---|
[9eefda] | 617 | } while (1); // (2)
|
---|
[cee0b57] | 618 |
|
---|
| 619 | // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
|
---|
[9eefda] | 620 | if ((Walker == DFS.Root) && (Binder == NULL))
|
---|
[cee0b57] | 621 | break;
|
---|
| 622 |
|
---|
[e138de] | 623 | DepthFirstSearchAnalysis_CheckForaNewComponent(this, Walker, DFS, LeafWalker);
|
---|
[174e0e] | 624 |
|
---|
[e138de] | 625 | DepthFirstSearchAnalysis_CleanRootStackDownTillWalker(this, Walker, Binder, DFS, LeafWalker);
|
---|
[174e0e] | 626 |
|
---|
[9eefda] | 627 | } while ((DFS.BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
|
---|
[cee0b57] | 628 |
|
---|
| 629 | // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
|
---|
[a67d19] | 630 | DoLog(0) && (Log() << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << DFS.CurrentGraphNr << "." << endl);
|
---|
[986ed3] | 631 | LeafWalker->Leaf->Output((ofstream *)&(Log() << Verbose(0)));
|
---|
[a67d19] | 632 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
[cee0b57] | 633 |
|
---|
| 634 | // step on to next root
|
---|
[9879f6] | 635 | while ((iter != end()) && ((*iter)->GraphNr != -1)) {
|
---|
| 636 | //Log() << Verbose(1) << "Current next subgraph root candidate is " << (*iter)->Name << "." << endl;
|
---|
| 637 | if ((*iter)->GraphNr != -1) // if already discovered, step on
|
---|
| 638 | iter++;
|
---|
[cee0b57] | 639 | }
|
---|
| 640 | }
|
---|
| 641 | // set cyclic bond criterium on "same LP" basis
|
---|
[266237] | 642 | CyclicBondAnalysis();
|
---|
| 643 |
|
---|
[e138de] | 644 | OutputGraphInfoPerAtom();
|
---|
[266237] | 645 |
|
---|
[e138de] | 646 | OutputGraphInfoPerBond();
|
---|
[266237] | 647 |
|
---|
| 648 | // free all and exit
|
---|
[e138de] | 649 | DepthFirstSearchAnalysis_Finalize(DFS);
|
---|
[a67d19] | 650 | DoLog(0) && (Log() << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl);
|
---|
[266237] | 651 | return SubGraphs;
|
---|
[9eefda] | 652 | }
|
---|
| 653 | ;
|
---|
[266237] | 654 |
|
---|
| 655 | /** Scans through all bonds and set bond::Cyclic to true where atom::LowpointNr of both ends is equal: LP criterion.
|
---|
| 656 | */
|
---|
[fa649a] | 657 | void molecule::CyclicBondAnalysis() const
|
---|
[266237] | 658 | {
|
---|
| 659 | NoCyclicBonds = 0;
|
---|
[9d83b6] | 660 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 661 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 662 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 663 | BondRunner != ListOfBonds.end();
|
---|
| 664 | ++BondRunner)
|
---|
[e08c46] | 665 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
| 666 | if ((*BondRunner)->rightatom->LowpointNr == (*BondRunner)->leftatom->LowpointNr) { // cyclic ??
|
---|
| 667 | (*BondRunner)->Cyclic = true;
|
---|
| 668 | NoCyclicBonds++;
|
---|
| 669 | }
|
---|
[9d83b6] | 670 | }
|
---|
[9eefda] | 671 | }
|
---|
| 672 | ;
|
---|
[cee0b57] | 673 |
|
---|
[266237] | 674 | /** Output graph information per atom.
|
---|
| 675 | * \param *out output stream
|
---|
| 676 | */
|
---|
[e138de] | 677 | void molecule::OutputGraphInfoPerAtom() const
|
---|
[266237] | 678 | {
|
---|
[a67d19] | 679 | DoLog(1) && (Log() << Verbose(1) << "Final graph info for each atom is:" << endl);
|
---|
[c743f8] | 680 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputGraphInfo));
|
---|
[9eefda] | 681 | }
|
---|
| 682 | ;
|
---|
[cee0b57] | 683 |
|
---|
[266237] | 684 | /** Output graph information per bond.
|
---|
| 685 | * \param *out output stream
|
---|
| 686 | */
|
---|
[e138de] | 687 | void molecule::OutputGraphInfoPerBond() const
|
---|
[266237] | 688 | {
|
---|
[a67d19] | 689 | DoLog(1) && (Log() << Verbose(1) << "Final graph info for each bond is:" << endl);
|
---|
[9d83b6] | 690 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 691 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 692 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 693 | BondRunner != ListOfBonds.end();
|
---|
| 694 | ++BondRunner)
|
---|
[e08c46] | 695 | if ((*BondRunner)->leftatom == *AtomRunner) {
|
---|
[9d83b6] | 696 | const bond *Binder = *BondRunner;
|
---|
[f9183b] | 697 | if (DoLog(2)) {
|
---|
| 698 | ostream &out = (Log() << Verbose(2));
|
---|
| 699 | out << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
|
---|
| 700 | out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
|
---|
| 701 | Binder->leftatom->OutputComponentNumber(&out);
|
---|
| 702 | out << " === ";
|
---|
| 703 | out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
|
---|
| 704 | Binder->rightatom->OutputComponentNumber(&out);
|
---|
| 705 | out << ">." << endl;
|
---|
| 706 | }
|
---|
[e08c46] | 707 | if (Binder->Cyclic) // cyclic ??
|
---|
| 708 | DoLog(3) && (Log() << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl);
|
---|
| 709 | }
|
---|
[9d83b6] | 710 | }
|
---|
[9eefda] | 711 | }
|
---|
| 712 | ;
|
---|
| 713 |
|
---|
| 714 | /** Initialise each vertex as white with no predecessor, empty queue, color Root lightgray.
|
---|
| 715 | * \param *out output stream for debugging
|
---|
| 716 | * \param &BFS accounting structure
|
---|
| 717 | * \param AtomCount number of entries in the array to allocate
|
---|
| 718 | */
|
---|
[e138de] | 719 | void InitializeBFSAccounting(struct BFSAccounting &BFS, int AtomCount)
|
---|
[9eefda] | 720 | {
|
---|
| 721 | BFS.AtomCount = AtomCount;
|
---|
[920c70] | 722 | BFS.PredecessorList = new atom*[AtomCount];
|
---|
| 723 | BFS.ShortestPathList = new int[AtomCount];
|
---|
| 724 | BFS.ColorList = new enum Shading[AtomCount];
|
---|
[a564be] | 725 | BFS.BFSStack = new std::deque<atom *> (AtomCount);
|
---|
| 726 | BFS.TouchedStack = new std::deque<atom *> (AtomCount);
|
---|
[9eefda] | 727 |
|
---|
[920c70] | 728 | for (int i = AtomCount; i--;) {
|
---|
[9eefda] | 729 | BFS.ShortestPathList[i] = -1;
|
---|
[920c70] | 730 | BFS.PredecessorList[i] = 0;
|
---|
[c27778] | 731 | BFS.ColorList[i] = white;
|
---|
[920c70] | 732 | }
|
---|
[cee0b57] | 733 | };
|
---|
| 734 |
|
---|
[9eefda] | 735 | /** Free's accounting structure.
|
---|
| 736 | * \param *out output stream for debugging
|
---|
| 737 | * \param &BFS accounting structure
|
---|
| 738 | */
|
---|
[e138de] | 739 | void FinalizeBFSAccounting(struct BFSAccounting &BFS)
|
---|
[9eefda] | 740 | {
|
---|
[920c70] | 741 | delete[](BFS.PredecessorList);
|
---|
| 742 | delete[](BFS.ShortestPathList);
|
---|
| 743 | delete[](BFS.ColorList);
|
---|
[9eefda] | 744 | delete (BFS.BFSStack);
|
---|
[c27778] | 745 | delete (BFS.TouchedStack);
|
---|
[9eefda] | 746 | BFS.AtomCount = 0;
|
---|
| 747 | };
|
---|
| 748 |
|
---|
| 749 | /** Clean the accounting structure.
|
---|
| 750 | * \param *out output stream for debugging
|
---|
| 751 | * \param &BFS accounting structure
|
---|
[ef9aae] | 752 | */
|
---|
[e138de] | 753 | void CleanBFSAccounting(struct BFSAccounting &BFS)
|
---|
[ef9aae] | 754 | {
|
---|
[9eefda] | 755 | atom *Walker = NULL;
|
---|
[a564be] | 756 | while (!BFS.TouchedStack->empty()) {
|
---|
| 757 | Walker = BFS.TouchedStack->front();
|
---|
| 758 | BFS.TouchedStack->pop_front();
|
---|
[9eefda] | 759 | BFS.PredecessorList[Walker->nr] = NULL;
|
---|
| 760 | BFS.ShortestPathList[Walker->nr] = -1;
|
---|
| 761 | BFS.ColorList[Walker->nr] = white;
|
---|
[ef9aae] | 762 | }
|
---|
| 763 | };
|
---|
| 764 |
|
---|
[9eefda] | 765 | /** Resets shortest path list and BFSStack.
|
---|
| 766 | * \param *out output stream for debugging
|
---|
| 767 | * \param *&Walker current node, pushed onto BFSAccounting::BFSStack and BFSAccounting::TouchedStack
|
---|
| 768 | * \param &BFS accounting structure
|
---|
| 769 | */
|
---|
[e138de] | 770 | void ResetBFSAccounting(atom *&Walker, struct BFSAccounting &BFS)
|
---|
[ef9aae] | 771 | {
|
---|
[9eefda] | 772 | BFS.ShortestPathList[Walker->nr] = 0;
|
---|
[a564be] | 773 | BFS.BFSStack->clear(); // start with empty BFS stack
|
---|
| 774 | BFS.BFSStack->push_front(Walker);
|
---|
| 775 | BFS.TouchedStack->push_front(Walker);
|
---|
[ef9aae] | 776 | };
|
---|
| 777 |
|
---|
[9eefda] | 778 | /** Performs a BFS from \a *Root, trying to find the same node and hence a cycle.
|
---|
| 779 | * \param *out output stream for debugging
|
---|
| 780 | * \param *&BackEdge the edge from root that we don't want to move along
|
---|
| 781 | * \param &BFS accounting structure
|
---|
| 782 | */
|
---|
[e138de] | 783 | void CyclicStructureAnalysis_CyclicBFSFromRootToRoot(bond *&BackEdge, struct BFSAccounting &BFS)
|
---|
[ef9aae] | 784 | {
|
---|
| 785 | atom *Walker = NULL;
|
---|
| 786 | atom *OtherAtom = NULL;
|
---|
[9eefda] | 787 | do { // look for Root
|
---|
[a564be] | 788 | ASSERT(!BFS.BFSStack->empty(), "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - BFS.BFSStack is empty!");
|
---|
| 789 | Walker = BFS.BFSStack->front();
|
---|
| 790 | BFS.BFSStack->pop_front();
|
---|
[a67d19] | 791 | DoLog(2) && (Log() << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *BFS.Root << "." << endl);
|
---|
[9d83b6] | 792 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 793 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 794 | Runner != ListOfBonds.end();
|
---|
| 795 | ++Runner) {
|
---|
[ef9aae] | 796 | if ((*Runner) != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
|
---|
| 797 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
[9eefda] | 798 | #ifdef ADDHYDROGEN
|
---|
[83f176] | 799 | if (OtherAtom->getType()->getAtomicNumber() != 1) {
|
---|
[9eefda] | 800 | #endif
|
---|
[68f03d] | 801 | DoLog(2) && (Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "." << endl);
|
---|
[9eefda] | 802 | if (BFS.ColorList[OtherAtom->nr] == white) {
|
---|
[a564be] | 803 | BFS.TouchedStack->push_front(OtherAtom);
|
---|
[9eefda] | 804 | BFS.ColorList[OtherAtom->nr] = lightgray;
|
---|
| 805 | BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
|
---|
| 806 | BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
|
---|
[68f03d] | 807 | DoLog(2) && (Log() << Verbose(2) << "Coloring OtherAtom " << OtherAtom->getName() << " lightgray, its predecessor is " << Walker->getName() << " and its Shortest Path is " << BFS.ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl);
|
---|
[9eefda] | 808 | //if (BFS.ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
|
---|
[a67d19] | 809 | DoLog(3) && (Log() << Verbose(3) << "Putting OtherAtom into queue." << endl);
|
---|
[a564be] | 810 | BFS.BFSStack->push_front(OtherAtom);
|
---|
[9eefda] | 811 | //}
|
---|
[ef9aae] | 812 | } else {
|
---|
[a67d19] | 813 | DoLog(3) && (Log() << Verbose(3) << "Not Adding, has already been visited." << endl);
|
---|
[ef9aae] | 814 | }
|
---|
[9eefda] | 815 | if (OtherAtom == BFS.Root)
|
---|
| 816 | break;
|
---|
| 817 | #ifdef ADDHYDROGEN
|
---|
| 818 | } else {
|
---|
[a67d19] | 819 | DoLog(2) && (Log() << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl);
|
---|
[9eefda] | 820 | BFS.ColorList[OtherAtom->nr] = black;
|
---|
| 821 | }
|
---|
| 822 | #endif
|
---|
[ef9aae] | 823 | } else {
|
---|
[a67d19] | 824 | DoLog(2) && (Log() << Verbose(2) << "Bond " << *(*Runner) << " not Visiting, is the back edge." << endl);
|
---|
[ef9aae] | 825 | }
|
---|
| 826 | }
|
---|
[9eefda] | 827 | BFS.ColorList[Walker->nr] = black;
|
---|
[68f03d] | 828 | DoLog(1) && (Log() << Verbose(1) << "Coloring Walker " << Walker->getName() << " black." << endl);
|
---|
[9eefda] | 829 | if (OtherAtom == BFS.Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
|
---|
[ef9aae] | 830 | // step through predecessor list
|
---|
| 831 | while (OtherAtom != BackEdge->rightatom) {
|
---|
[9eefda] | 832 | if (!OtherAtom->GetTrueFather()->IsCyclic) // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
|
---|
[ef9aae] | 833 | break;
|
---|
| 834 | else
|
---|
[9eefda] | 835 | OtherAtom = BFS.PredecessorList[OtherAtom->nr];
|
---|
[ef9aae] | 836 | }
|
---|
| 837 | if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
|
---|
[a67d19] | 838 | DoLog(3) && (Log() << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl);
|
---|
[ef9aae] | 839 | do {
|
---|
[a564be] | 840 | ASSERT(!BFS.TouchedStack->empty(), "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - BFS.TouchedStack is empty!");
|
---|
| 841 | OtherAtom = BFS.TouchedStack->front();
|
---|
| 842 | BFS.TouchedStack->pop_front();
|
---|
[9eefda] | 843 | if (BFS.PredecessorList[OtherAtom->nr] == Walker) {
|
---|
[a67d19] | 844 | DoLog(4) && (Log() << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl);
|
---|
[9eefda] | 845 | BFS.PredecessorList[OtherAtom->nr] = NULL;
|
---|
| 846 | BFS.ShortestPathList[OtherAtom->nr] = -1;
|
---|
| 847 | BFS.ColorList[OtherAtom->nr] = white;
|
---|
[a564be] | 848 | // rats ... deque has no find()
|
---|
| 849 | std::deque<atom *>::iterator iter = find(
|
---|
| 850 | BFS.BFSStack->begin(),
|
---|
| 851 | BFS.BFSStack->end(),
|
---|
| 852 | OtherAtom);
|
---|
| 853 | ASSERT(iter != BFS.BFSStack->end(),
|
---|
| 854 | "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - can't find "+toString(*OtherAtom)+" on stack!");
|
---|
| 855 | BFS.BFSStack->erase(iter);
|
---|
[ef9aae] | 856 | }
|
---|
[a564be] | 857 | } while ((!BFS.TouchedStack->empty()) && (BFS.PredecessorList[OtherAtom->nr] == NULL));
|
---|
| 858 | BFS.TouchedStack->push_front(OtherAtom); // last was wrongly popped
|
---|
[ef9aae] | 859 | OtherAtom = BackEdge->rightatom; // set to not Root
|
---|
| 860 | } else
|
---|
[9eefda] | 861 | OtherAtom = BFS.Root;
|
---|
[ef9aae] | 862 | }
|
---|
[a564be] | 863 | } while ((!BFS.BFSStack->empty()) && (OtherAtom != BFS.Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
|
---|
[ef9aae] | 864 | };
|
---|
| 865 |
|
---|
[9eefda] | 866 | /** Climb back the BFSAccounting::PredecessorList and find cycle members.
|
---|
| 867 | * \param *out output stream for debugging
|
---|
| 868 | * \param *&OtherAtom
|
---|
| 869 | * \param *&BackEdge denotes the edge we did not want to travel along when doing CyclicBFSFromRootToRoot()
|
---|
| 870 | * \param &BFS accounting structure
|
---|
| 871 | * \param *&MinimumRingSize minimum distance from this node possible without encountering oneself, set on return for each atom
|
---|
| 872 | * \param &MinRingSize global minimum distance from one node without encountering oneself, set on return
|
---|
| 873 | */
|
---|
[e138de] | 874 | void CyclicStructureAnalysis_RetrieveCycleMembers(atom *&OtherAtom, bond *&BackEdge, struct BFSAccounting &BFS, int *&MinimumRingSize, int &MinRingSize)
|
---|
[ef9aae] | 875 | {
|
---|
| 876 | atom *Walker = NULL;
|
---|
| 877 | int NumCycles = 0;
|
---|
| 878 | int RingSize = -1;
|
---|
| 879 |
|
---|
[9eefda] | 880 | if (OtherAtom == BFS.Root) {
|
---|
[ef9aae] | 881 | // now climb back the predecessor list and thus find the cycle members
|
---|
| 882 | NumCycles++;
|
---|
| 883 | RingSize = 1;
|
---|
[9eefda] | 884 | BFS.Root->GetTrueFather()->IsCyclic = true;
|
---|
[a67d19] | 885 | DoLog(1) && (Log() << Verbose(1) << "Found ring contains: ");
|
---|
[9eefda] | 886 | Walker = BFS.Root;
|
---|
[ef9aae] | 887 | while (Walker != BackEdge->rightatom) {
|
---|
[68f03d] | 888 | DoLog(0) && (Log() << Verbose(0) << Walker->getName() << " <-> ");
|
---|
[9eefda] | 889 | Walker = BFS.PredecessorList[Walker->nr];
|
---|
[ef9aae] | 890 | Walker->GetTrueFather()->IsCyclic = true;
|
---|
| 891 | RingSize++;
|
---|
| 892 | }
|
---|
[68f03d] | 893 | DoLog(0) && (Log() << Verbose(0) << Walker->getName() << " with a length of " << RingSize << "." << endl << endl);
|
---|
[ef9aae] | 894 | // walk through all and set MinimumRingSize
|
---|
[9eefda] | 895 | Walker = BFS.Root;
|
---|
[ef9aae] | 896 | MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
|
---|
| 897 | while (Walker != BackEdge->rightatom) {
|
---|
[9eefda] | 898 | Walker = BFS.PredecessorList[Walker->nr];
|
---|
[ef9aae] | 899 | if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
|
---|
| 900 | MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
|
---|
| 901 | }
|
---|
| 902 | if ((RingSize < MinRingSize) || (MinRingSize == -1))
|
---|
| 903 | MinRingSize = RingSize;
|
---|
| 904 | } else {
|
---|
[c27778] | 905 | DoLog(1) && (Log() << Verbose(1) << "No ring containing " << *BFS.Root << " with length equal to or smaller than " << MinimumRingSize[BFS.Root->GetTrueFather()->nr] << " found." << endl);
|
---|
[ef9aae] | 906 | }
|
---|
| 907 | };
|
---|
| 908 |
|
---|
[9eefda] | 909 | /** From a given node performs a BFS to touch the next cycle, for whose nodes \a *&MinimumRingSize is set and set it accordingly.
|
---|
| 910 | * \param *out output stream for debugging
|
---|
| 911 | * \param *&Root node to look for closest cycle from, i.e. \a *&MinimumRingSize is set for this node
|
---|
| 912 | * \param *&MinimumRingSize minimum distance from this node possible without encountering oneself, set on return for each atom
|
---|
| 913 | * \param AtomCount number of nodes in graph
|
---|
| 914 | */
|
---|
[e138de] | 915 | void CyclicStructureAnalysis_BFSToNextCycle(atom *&Root, atom *&Walker, int *&MinimumRingSize, int AtomCount)
|
---|
[ef9aae] | 916 | {
|
---|
[9eefda] | 917 | struct BFSAccounting BFS;
|
---|
[ef9aae] | 918 | atom *OtherAtom = Walker;
|
---|
| 919 |
|
---|
[e138de] | 920 | InitializeBFSAccounting(BFS, AtomCount);
|
---|
[ef9aae] | 921 |
|
---|
[e138de] | 922 | ResetBFSAccounting(Walker, BFS);
|
---|
[9eefda] | 923 | while (OtherAtom != NULL) { // look for Root
|
---|
[a564be] | 924 | ASSERT(!BFS.BFSStack->empty(), "CyclicStructureAnalysis_BFSToNextCycle() - BFS.BFSStack is empty!");
|
---|
| 925 | Walker = BFS.BFSStack->front();
|
---|
| 926 | BFS.BFSStack->pop_front();
|
---|
[e138de] | 927 | //Log() << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
|
---|
[9d83b6] | 928 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 929 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 930 | Runner != ListOfBonds.end();
|
---|
| 931 | ++Runner) {
|
---|
[9eefda] | 932 | // "removed (*Runner) != BackEdge) || " from next if, is u
|
---|
[9d83b6] | 933 | if ((ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
|
---|
[ef9aae] | 934 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
[e138de] | 935 | //Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
|
---|
[9eefda] | 936 | if (BFS.ColorList[OtherAtom->nr] == white) {
|
---|
[a564be] | 937 | BFS.TouchedStack->push_front(OtherAtom);
|
---|
[9eefda] | 938 | BFS.ColorList[OtherAtom->nr] = lightgray;
|
---|
| 939 | BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
|
---|
| 940 | BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
|
---|
[e138de] | 941 | //Log() << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
|
---|
[ef9aae] | 942 | if (OtherAtom->GetTrueFather()->IsCyclic) { // if the other atom is connected to a ring
|
---|
[9eefda] | 943 | MinimumRingSize[Root->GetTrueFather()->nr] = BFS.ShortestPathList[OtherAtom->nr] + MinimumRingSize[OtherAtom->GetTrueFather()->nr];
|
---|
[ef9aae] | 944 | OtherAtom = NULL; //break;
|
---|
| 945 | break;
|
---|
| 946 | } else
|
---|
[a564be] | 947 | BFS.BFSStack->push_front(OtherAtom);
|
---|
[ef9aae] | 948 | } else {
|
---|
[e138de] | 949 | //Log() << Verbose(3) << "Not Adding, has already been visited." << endl;
|
---|
[ef9aae] | 950 | }
|
---|
| 951 | } else {
|
---|
[e138de] | 952 | //Log() << Verbose(3) << "Not Visiting, is a back edge." << endl;
|
---|
[ef9aae] | 953 | }
|
---|
| 954 | }
|
---|
[9eefda] | 955 | BFS.ColorList[Walker->nr] = black;
|
---|
[e138de] | 956 | //Log() << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
|
---|
[ef9aae] | 957 | }
|
---|
| 958 | //CleanAccountingLists(TouchedStack, PredecessorList, ShortestPathList, ColorList);
|
---|
| 959 |
|
---|
[e138de] | 960 | FinalizeBFSAccounting(BFS);
|
---|
[9eefda] | 961 | }
|
---|
| 962 | ;
|
---|
[ef9aae] | 963 |
|
---|
[9eefda] | 964 | /** All nodes that are not in cycles get assigned a \a *&MinimumRingSizeby BFS to next cycle.
|
---|
| 965 | * \param *out output stream for debugging
|
---|
| 966 | * \param *&MinimumRingSize array with minimum distance without encountering onself for each atom
|
---|
| 967 | * \param &MinRingSize global minium distance
|
---|
| 968 | * \param &NumCyles number of cycles in graph
|
---|
| 969 | * \param *mol molecule with atoms
|
---|
| 970 | */
|
---|
[e138de] | 971 | void CyclicStructureAnalysis_AssignRingSizetoNonCycleMembers(int *&MinimumRingSize, int &MinRingSize, int &NumCycles, const molecule * const mol)
|
---|
[ef9aae] | 972 | {
|
---|
[9eefda] | 973 | atom *Root = NULL;
|
---|
[ef9aae] | 974 | atom *Walker = NULL;
|
---|
| 975 | if (MinRingSize != -1) { // if rings are present
|
---|
| 976 | // go over all atoms
|
---|
[9879f6] | 977 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
| 978 | Root = *iter;
|
---|
[ef9aae] | 979 |
|
---|
[ea7176] | 980 | if (MinimumRingSize[Root->GetTrueFather()->nr] == mol->getAtomCount()) { // check whether MinimumRingSize is set, if not BFS to next where it is
|
---|
[ef9aae] | 981 | Walker = Root;
|
---|
| 982 |
|
---|
[e138de] | 983 | //Log() << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
|
---|
[ea7176] | 984 | CyclicStructureAnalysis_BFSToNextCycle(Root, Walker, MinimumRingSize, mol->getAtomCount());
|
---|
[ef9aae] | 985 |
|
---|
| 986 | }
|
---|
[a67d19] | 987 | DoLog(1) && (Log() << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl);
|
---|
[ef9aae] | 988 | }
|
---|
[a67d19] | 989 | DoLog(1) && (Log() << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl);
|
---|
[ef9aae] | 990 | } else
|
---|
[a67d19] | 991 | DoLog(1) && (Log() << Verbose(1) << "No rings were detected in the molecular structure." << endl);
|
---|
[9eefda] | 992 | }
|
---|
| 993 | ;
|
---|
[ef9aae] | 994 |
|
---|
[cee0b57] | 995 | /** Analyses the cycles found and returns minimum of all cycle lengths.
|
---|
| 996 | * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root,
|
---|
| 997 | * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
|
---|
| 998 | * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
|
---|
| 999 | * as cyclic and print out the cycles.
|
---|
| 1000 | * \param *out output stream for debugging
|
---|
| 1001 | * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
|
---|
| 1002 | * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
|
---|
| 1003 | * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
|
---|
| 1004 | */
|
---|
[a564be] | 1005 | void molecule::CyclicStructureAnalysis(std::deque<bond *> * BackEdgeStack, int *&MinimumRingSize) const
|
---|
[cee0b57] | 1006 | {
|
---|
[9eefda] | 1007 | struct BFSAccounting BFS;
|
---|
[ef9aae] | 1008 | atom *Walker = NULL;
|
---|
| 1009 | atom *OtherAtom = NULL;
|
---|
| 1010 | bond *BackEdge = NULL;
|
---|
| 1011 | int NumCycles = 0;
|
---|
| 1012 | int MinRingSize = -1;
|
---|
[cee0b57] | 1013 |
|
---|
[ea7176] | 1014 | InitializeBFSAccounting(BFS, getAtomCount());
|
---|
[cee0b57] | 1015 |
|
---|
[e138de] | 1016 | //Log() << Verbose(1) << "Back edge list - ";
|
---|
[99593f] | 1017 | //BackEdgeStack->Output(out);
|
---|
[cee0b57] | 1018 |
|
---|
[a67d19] | 1019 | DoLog(1) && (Log() << Verbose(1) << "Analysing cycles ... " << endl);
|
---|
[cee0b57] | 1020 | NumCycles = 0;
|
---|
[a564be] | 1021 | while (!BackEdgeStack->empty()) {
|
---|
| 1022 | BackEdge = BackEdgeStack->front();
|
---|
| 1023 | BackEdgeStack->pop_front();
|
---|
[cee0b57] | 1024 | // this is the target
|
---|
[9eefda] | 1025 | BFS.Root = BackEdge->leftatom;
|
---|
[cee0b57] | 1026 | // this is the source point
|
---|
| 1027 | Walker = BackEdge->rightatom;
|
---|
| 1028 |
|
---|
[e138de] | 1029 | ResetBFSAccounting(Walker, BFS);
|
---|
[cee0b57] | 1030 |
|
---|
[a67d19] | 1031 | DoLog(1) && (Log() << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl);
|
---|
[ef9aae] | 1032 | OtherAtom = NULL;
|
---|
[e138de] | 1033 | CyclicStructureAnalysis_CyclicBFSFromRootToRoot(BackEdge, BFS);
|
---|
[cee0b57] | 1034 |
|
---|
[e138de] | 1035 | CyclicStructureAnalysis_RetrieveCycleMembers(OtherAtom, BackEdge, BFS, MinimumRingSize, MinRingSize);
|
---|
[cee0b57] | 1036 |
|
---|
[e138de] | 1037 | CleanBFSAccounting(BFS);
|
---|
[ef9aae] | 1038 | }
|
---|
[e138de] | 1039 | FinalizeBFSAccounting(BFS);
|
---|
[ef9aae] | 1040 |
|
---|
[e138de] | 1041 | CyclicStructureAnalysis_AssignRingSizetoNonCycleMembers(MinimumRingSize, MinRingSize, NumCycles, this);
|
---|
[fa649a] | 1042 | };
|
---|
[cee0b57] | 1043 |
|
---|
| 1044 | /** Sets the next component number.
|
---|
| 1045 | * This is O(N) as the number of bonds per atom is bound.
|
---|
| 1046 | * \param *vertex atom whose next atom::*ComponentNr is to be set
|
---|
| 1047 | * \param nr number to use
|
---|
| 1048 | */
|
---|
[fa649a] | 1049 | void molecule::SetNextComponentNumber(atom *vertex, int nr) const
|
---|
[cee0b57] | 1050 | {
|
---|
[9eefda] | 1051 | size_t i = 0;
|
---|
[cee0b57] | 1052 | if (vertex != NULL) {
|
---|
[9d83b6] | 1053 | const BondList& ListOfBonds = vertex->getListOfBonds();
|
---|
| 1054 | for (; i < ListOfBonds.size(); i++) {
|
---|
[9eefda] | 1055 | if (vertex->ComponentNr[i] == -1) { // check if not yet used
|
---|
[cee0b57] | 1056 | vertex->ComponentNr[i] = nr;
|
---|
| 1057 | break;
|
---|
[9eefda] | 1058 | } else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
|
---|
| 1059 | break; // breaking here will not cause error!
|
---|
[cee0b57] | 1060 | }
|
---|
[9d83b6] | 1061 | if (i == ListOfBonds.size()) {
|
---|
[58ed4a] | 1062 | DoeLog(0) && (eLog()<< Verbose(0) << "Error: All Component entries are already occupied!" << endl);
|
---|
[e359a8] | 1063 | performCriticalExit();
|
---|
| 1064 | }
|
---|
| 1065 | } else {
|
---|
[58ed4a] | 1066 | DoeLog(0) && (eLog()<< Verbose(0) << "Error: Given vertex is NULL!" << endl);
|
---|
[e359a8] | 1067 | performCriticalExit();
|
---|
| 1068 | }
|
---|
[9eefda] | 1069 | }
|
---|
| 1070 | ;
|
---|
[cee0b57] | 1071 |
|
---|
| 1072 | /** Returns next unused bond for this atom \a *vertex or NULL of none exists.
|
---|
| 1073 | * \param *vertex atom to regard
|
---|
| 1074 | * \return bond class or NULL
|
---|
| 1075 | */
|
---|
[fa649a] | 1076 | bond * molecule::FindNextUnused(atom *vertex) const
|
---|
[cee0b57] | 1077 | {
|
---|
[9d83b6] | 1078 | const BondList& ListOfBonds = vertex->getListOfBonds();
|
---|
| 1079 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 1080 | Runner != ListOfBonds.end();
|
---|
| 1081 | ++Runner)
|
---|
[266237] | 1082 | if ((*Runner)->IsUsed() == white)
|
---|
[9eefda] | 1083 | return ((*Runner));
|
---|
[cee0b57] | 1084 | return NULL;
|
---|
[9eefda] | 1085 | }
|
---|
| 1086 | ;
|
---|
[cee0b57] | 1087 |
|
---|
| 1088 | /** Resets bond::Used flag of all bonds in this molecule.
|
---|
| 1089 | * \return true - success, false - -failure
|
---|
| 1090 | */
|
---|
[fa649a] | 1091 | void molecule::ResetAllBondsToUnused() const
|
---|
[cee0b57] | 1092 | {
|
---|
[9d83b6] | 1093 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
|
---|
| 1094 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 1095 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 1096 | BondRunner != ListOfBonds.end();
|
---|
| 1097 | ++BondRunner)
|
---|
[e08c46] | 1098 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
| 1099 | (*BondRunner)->ResetUsed();
|
---|
[9d83b6] | 1100 | }
|
---|
[9eefda] | 1101 | }
|
---|
| 1102 | ;
|
---|
[cee0b57] | 1103 |
|
---|
| 1104 | /** Output a list of flags, stating whether the bond was visited or not.
|
---|
| 1105 | * \param *out output stream for debugging
|
---|
| 1106 | * \param *list
|
---|
| 1107 | */
|
---|
[e138de] | 1108 | void OutputAlreadyVisited(int *list)
|
---|
[cee0b57] | 1109 | {
|
---|
[a67d19] | 1110 | DoLog(4) && (Log() << Verbose(4) << "Already Visited Bonds:\t");
|
---|
[9eefda] | 1111 | for (int i = 1; i <= list[0]; i++)
|
---|
[a67d19] | 1112 | DoLog(0) && (Log() << Verbose(0) << list[i] << " ");
|
---|
| 1113 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
[9eefda] | 1114 | }
|
---|
| 1115 | ;
|
---|
[cee0b57] | 1116 |
|
---|
| 1117 | /** Storing the bond structure of a molecule to file.
|
---|
| 1118 | * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
|
---|
[35b698] | 1119 | * \param &filename name of file
|
---|
| 1120 | * \param path path to file, defaults to empty
|
---|
[cee0b57] | 1121 | * \return true - file written successfully, false - writing failed
|
---|
| 1122 | */
|
---|
[e4afb4] | 1123 | bool molecule::StoreAdjacencyToFile(std::string filename, std::string path)
|
---|
[cee0b57] | 1124 | {
|
---|
| 1125 | ofstream AdjacencyFile;
|
---|
[35b698] | 1126 | string line;
|
---|
[cee0b57] | 1127 | bool status = true;
|
---|
| 1128 |
|
---|
[35b698] | 1129 | if (path != "")
|
---|
| 1130 | line = path + "/" + filename;
|
---|
[8ab0407] | 1131 | else
|
---|
[35b698] | 1132 | line = filename;
|
---|
| 1133 | AdjacencyFile.open(line.c_str(), ios::out);
|
---|
[acf800] | 1134 | DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
|
---|
[35b698] | 1135 | if (AdjacencyFile.good()) {
|
---|
[1f1b23] | 1136 | AdjacencyFile << "m\tn" << endl;
|
---|
[00ef5c] | 1137 | for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputAdjacency),&AdjacencyFile));
|
---|
[cee0b57] | 1138 | AdjacencyFile.close();
|
---|
[acf800] | 1139 | DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
|
---|
[cee0b57] | 1140 | } else {
|
---|
[35b698] | 1141 | DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
|
---|
[cee0b57] | 1142 | status = false;
|
---|
| 1143 | }
|
---|
| 1144 |
|
---|
| 1145 | return status;
|
---|
[9eefda] | 1146 | }
|
---|
| 1147 | ;
|
---|
[cee0b57] | 1148 |
|
---|
[1f1b23] | 1149 | /** Storing the bond structure of a molecule to file.
|
---|
| 1150 | * Simply stores Atom::nr and then the Atom::nr of all bond partners, one per line.
|
---|
[35b698] | 1151 | * \param &filename name of file
|
---|
| 1152 | * \param path path to file, defaults to empty
|
---|
[1f1b23] | 1153 | * \return true - file written successfully, false - writing failed
|
---|
| 1154 | */
|
---|
[e4afb4] | 1155 | bool molecule::StoreBondsToFile(std::string filename, std::string path)
|
---|
[1f1b23] | 1156 | {
|
---|
| 1157 | ofstream BondFile;
|
---|
[35b698] | 1158 | string line;
|
---|
[1f1b23] | 1159 | bool status = true;
|
---|
| 1160 |
|
---|
[35b698] | 1161 | if (path != "")
|
---|
| 1162 | line = path + "/" + filename;
|
---|
[8ab0407] | 1163 | else
|
---|
[35b698] | 1164 | line = filename;
|
---|
| 1165 | BondFile.open(line.c_str(), ios::out);
|
---|
[acf800] | 1166 | DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
|
---|
[35b698] | 1167 | if (BondFile.good()) {
|
---|
[1f1b23] | 1168 | BondFile << "m\tn" << endl;
|
---|
[00ef5c] | 1169 | for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile));
|
---|
[1f1b23] | 1170 | BondFile.close();
|
---|
[acf800] | 1171 | DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
|
---|
[1f1b23] | 1172 | } else {
|
---|
[35b698] | 1173 | DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
|
---|
[1f1b23] | 1174 | status = false;
|
---|
| 1175 | }
|
---|
| 1176 |
|
---|
| 1177 | return status;
|
---|
| 1178 | }
|
---|
| 1179 | ;
|
---|
| 1180 |
|
---|
[35b698] | 1181 | bool CheckAdjacencyFileAgainstMolecule_Init(std::string &path, ifstream &File, int *&CurrentBonds)
|
---|
[ba4170] | 1182 | {
|
---|
[35b698] | 1183 | string filename;
|
---|
| 1184 | filename = path + ADJACENCYFILE;
|
---|
| 1185 | File.open(filename.c_str(), ios::out);
|
---|
[0de7e8] | 1186 | DoLog(1) && (Log() << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... " << endl);
|
---|
[35b698] | 1187 | if (File.fail())
|
---|
[ba4170] | 1188 | return false;
|
---|
| 1189 |
|
---|
| 1190 | // allocate storage structure
|
---|
[1d5afa5] | 1191 | CurrentBonds = new int[MAXBONDS]; // contains parsed bonds of current atom
|
---|
| 1192 | for(int i=0;i<MAXBONDS;i++)
|
---|
[920c70] | 1193 | CurrentBonds[i] = 0;
|
---|
[ba4170] | 1194 | return true;
|
---|
[9eefda] | 1195 | }
|
---|
| 1196 | ;
|
---|
[ba4170] | 1197 |
|
---|
[e138de] | 1198 | void CheckAdjacencyFileAgainstMolecule_Finalize(ifstream &File, int *&CurrentBonds)
|
---|
[ba4170] | 1199 | {
|
---|
| 1200 | File.close();
|
---|
| 1201 | File.clear();
|
---|
[920c70] | 1202 | delete[](CurrentBonds);
|
---|
[9eefda] | 1203 | }
|
---|
| 1204 | ;
|
---|
[ba4170] | 1205 |
|
---|
[e138de] | 1206 | void CheckAdjacencyFileAgainstMolecule_CompareBonds(bool &status, int &NonMatchNumber, atom *&Walker, size_t &CurrentBondsOfAtom, int AtomNr, int *&CurrentBonds, atom **ListOfAtoms)
|
---|
[ba4170] | 1207 | {
|
---|
| 1208 | size_t j = 0;
|
---|
| 1209 | int id = -1;
|
---|
| 1210 |
|
---|
[e138de] | 1211 | //Log() << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
|
---|
[9d83b6] | 1212 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 1213 | if (CurrentBondsOfAtom == ListOfBonds.size()) {
|
---|
| 1214 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 1215 | Runner != ListOfBonds.end();
|
---|
| 1216 | ++Runner) {
|
---|
[ba4170] | 1217 | id = (*Runner)->GetOtherAtom(Walker)->nr;
|
---|
| 1218 | j = 0;
|
---|
[9eefda] | 1219 | for (; (j < CurrentBondsOfAtom) && (CurrentBonds[j++] != id);)
|
---|
[ba4170] | 1220 | ; // check against all parsed bonds
|
---|
[9eefda] | 1221 | if (CurrentBonds[j - 1] != id) { // no match ? Then mark in ListOfAtoms
|
---|
[ba4170] | 1222 | ListOfAtoms[AtomNr] = NULL;
|
---|
| 1223 | NonMatchNumber++;
|
---|
| 1224 | status = false;
|
---|
[0de7e8] | 1225 | DoeLog(2) && (eLog() << Verbose(2) << id << " can not be found in list." << endl);
|
---|
[ba4170] | 1226 | } else {
|
---|
[0de7e8] | 1227 | //Log() << Verbose(0) << "[" << id << "]\t";
|
---|
[ba4170] | 1228 | }
|
---|
| 1229 | }
|
---|
[e138de] | 1230 | //Log() << Verbose(0) << endl;
|
---|
[ba4170] | 1231 | } else {
|
---|
[9d83b6] | 1232 | DoLog(0) && (Log() << Verbose(0) << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << ListOfBonds.size() << "." << endl);
|
---|
[ba4170] | 1233 | status = false;
|
---|
| 1234 | }
|
---|
[9eefda] | 1235 | }
|
---|
| 1236 | ;
|
---|
[ba4170] | 1237 |
|
---|
[cee0b57] | 1238 | /** Checks contents of adjacency file against bond structure in structure molecule.
|
---|
| 1239 | * \param *out output stream for debugging
|
---|
| 1240 | * \param *path path to file
|
---|
| 1241 | * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
|
---|
| 1242 | * \return true - structure is equal, false - not equivalence
|
---|
| 1243 | */
|
---|
[35b698] | 1244 | bool molecule::CheckAdjacencyFileAgainstMolecule(std::string &path, atom **ListOfAtoms)
|
---|
[cee0b57] | 1245 | {
|
---|
| 1246 | ifstream File;
|
---|
| 1247 | bool status = true;
|
---|
[266237] | 1248 | atom *Walker = NULL;
|
---|
[ba4170] | 1249 | int *CurrentBonds = NULL;
|
---|
[9eefda] | 1250 | int NonMatchNumber = 0; // will number of atoms with differing bond structure
|
---|
[ba4170] | 1251 | size_t CurrentBondsOfAtom = -1;
|
---|
[0de7e8] | 1252 | const int AtomCount = getAtomCount();
|
---|
[cee0b57] | 1253 |
|
---|
[e138de] | 1254 | if (!CheckAdjacencyFileAgainstMolecule_Init(path, File, CurrentBonds)) {
|
---|
[a67d19] | 1255 | DoLog(1) && (Log() << Verbose(1) << "Adjacency file not found." << endl);
|
---|
[ba4170] | 1256 | return true;
|
---|
| 1257 | }
|
---|
| 1258 |
|
---|
[920c70] | 1259 | char buffer[MAXSTRINGSIZE];
|
---|
[1d5afa5] | 1260 | int tmp;
|
---|
[ba4170] | 1261 | // Parse the file line by line and count the bonds
|
---|
| 1262 | while (!File.eof()) {
|
---|
| 1263 | File.getline(buffer, MAXSTRINGSIZE);
|
---|
| 1264 | stringstream line;
|
---|
| 1265 | line.str(buffer);
|
---|
| 1266 | int AtomNr = -1;
|
---|
| 1267 | line >> AtomNr;
|
---|
| 1268 | CurrentBondsOfAtom = -1; // we count one too far due to line end
|
---|
| 1269 | // parse into structure
|
---|
[0de7e8] | 1270 | if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
|
---|
[ba4170] | 1271 | Walker = ListOfAtoms[AtomNr];
|
---|
[1d5afa5] | 1272 | while (line >> ws >> tmp) {
|
---|
| 1273 | std::cout << "Recognized bond partner " << tmp << std::endl;
|
---|
| 1274 | CurrentBonds[++CurrentBondsOfAtom] = tmp;
|
---|
| 1275 | ASSERT(CurrentBondsOfAtom < MAXBONDS,
|
---|
| 1276 | "molecule::CheckAdjacencyFileAgainstMolecule() - encountered more bonds than allowed: "
|
---|
| 1277 | +toString(CurrentBondsOfAtom)+" >= "+toString(MAXBONDS)+"!");
|
---|
| 1278 | }
|
---|
[ba4170] | 1279 | // compare against present bonds
|
---|
[e138de] | 1280 | CheckAdjacencyFileAgainstMolecule_CompareBonds(status, NonMatchNumber, Walker, CurrentBondsOfAtom, AtomNr, CurrentBonds, ListOfAtoms);
|
---|
[0de7e8] | 1281 | } else {
|
---|
| 1282 | if (AtomNr != -1)
|
---|
| 1283 | DoeLog(2) && (eLog() << Verbose(2) << AtomNr << " is not valid in the range of ids [" << 0 << "," << AtomCount << ")." << endl);
|
---|
[ba4170] | 1284 | }
|
---|
[cee0b57] | 1285 | }
|
---|
[e138de] | 1286 | CheckAdjacencyFileAgainstMolecule_Finalize(File, CurrentBonds);
|
---|
[cee0b57] | 1287 |
|
---|
[ba4170] | 1288 | if (status) { // if equal we parse the KeySetFile
|
---|
[a67d19] | 1289 | DoLog(1) && (Log() << Verbose(1) << "done: Equal." << endl);
|
---|
[ba4170] | 1290 | } else
|
---|
[a67d19] | 1291 | DoLog(1) && (Log() << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl);
|
---|
[cee0b57] | 1292 | return status;
|
---|
[9eefda] | 1293 | }
|
---|
| 1294 | ;
|
---|
[cee0b57] | 1295 |
|
---|
| 1296 | /** Picks from a global stack with all back edges the ones in the fragment.
|
---|
| 1297 | * \param *out output stream for debugging
|
---|
| 1298 | * \param **ListOfLocalAtoms array of father atom::nr to local atom::nr (reverse of atom::father)
|
---|
| 1299 | * \param *ReferenceStack stack with all the back egdes
|
---|
| 1300 | * \param *LocalStack stack to be filled
|
---|
| 1301 | * \return true - everything ok, false - ReferenceStack was empty
|
---|
| 1302 | */
|
---|
[a564be] | 1303 | bool molecule::PickLocalBackEdges(atom **ListOfLocalAtoms, std::deque<bond *> *&ReferenceStack, std::deque<bond *> *&LocalStack) const
|
---|
[cee0b57] | 1304 | {
|
---|
| 1305 | bool status = true;
|
---|
[a564be] | 1306 | if (ReferenceStack->empty()) {
|
---|
[a67d19] | 1307 | DoLog(1) && (Log() << Verbose(1) << "ReferenceStack is empty!" << endl);
|
---|
[cee0b57] | 1308 | return false;
|
---|
| 1309 | }
|
---|
[a564be] | 1310 | bond *Binder = ReferenceStack->front();
|
---|
| 1311 | ReferenceStack->pop_front();
|
---|
[9eefda] | 1312 | bond *FirstBond = Binder; // mark the first bond, so that we don't loop through the stack indefinitely
|
---|
[cee0b57] | 1313 | atom *Walker = NULL, *OtherAtom = NULL;
|
---|
[a564be] | 1314 | ReferenceStack->push_front(Binder);
|
---|
[cee0b57] | 1315 |
|
---|
[9eefda] | 1316 | do { // go through all bonds and push local ones
|
---|
| 1317 | Walker = ListOfLocalAtoms[Binder->leftatom->nr]; // get one atom in the reference molecule
|
---|
[9d83b6] | 1318 | if (Walker != NULL) { // if this Walker exists in the subgraph ...
|
---|
| 1319 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 1320 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 1321 | Runner != ListOfBonds.end();
|
---|
| 1322 | ++Runner) {
|
---|
[266237] | 1323 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
| 1324 | if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->nr]) { // found the bond
|
---|
[a564be] | 1325 | LocalStack->push_front((*Runner));
|
---|
[a67d19] | 1326 | DoLog(3) && (Log() << Verbose(3) << "Found local edge " << *(*Runner) << "." << endl);
|
---|
[cee0b57] | 1327 | break;
|
---|
| 1328 | }
|
---|
| 1329 | }
|
---|
[9d83b6] | 1330 | }
|
---|
[a564be] | 1331 | ASSERT(!ReferenceStack->empty(), "molecule::PickLocalBackEdges() - ReferenceStack is empty!");
|
---|
| 1332 | Binder = ReferenceStack->front(); // loop the stack for next item
|
---|
| 1333 | ReferenceStack->pop_front();
|
---|
[a67d19] | 1334 | DoLog(3) && (Log() << Verbose(3) << "Current candidate edge " << Binder << "." << endl);
|
---|
[a564be] | 1335 | ReferenceStack->push_front(Binder);
|
---|
[cee0b57] | 1336 | } while (FirstBond != Binder);
|
---|
| 1337 |
|
---|
| 1338 | return status;
|
---|
[9eefda] | 1339 | }
|
---|
| 1340 | ;
|
---|
[ce7cc5] | 1341 |
|
---|
| 1342 | void BreadthFirstSearchAdd_Init(struct BFSAccounting &BFS, atom *&Root, int AtomCount, int BondOrder, atom **AddedAtomList = NULL)
|
---|
| 1343 | {
|
---|
| 1344 | BFS.AtomCount = AtomCount;
|
---|
| 1345 | BFS.BondOrder = BondOrder;
|
---|
[920c70] | 1346 | BFS.PredecessorList = new atom*[AtomCount];
|
---|
| 1347 | BFS.ShortestPathList = new int[AtomCount];
|
---|
| 1348 | BFS.ColorList = new enum Shading[AtomCount];
|
---|
[a564be] | 1349 | BFS.BFSStack = new std::deque<atom *> (AtomCount);
|
---|
[ce7cc5] | 1350 |
|
---|
| 1351 | BFS.Root = Root;
|
---|
[a564be] | 1352 | BFS.BFSStack->clear();
|
---|
| 1353 | BFS.BFSStack->push_front(Root);
|
---|
[ce7cc5] | 1354 |
|
---|
| 1355 | // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
|
---|
[9eefda] | 1356 | for (int i = AtomCount; i--;) {
|
---|
[920c70] | 1357 | BFS.PredecessorList[i] = NULL;
|
---|
[ce7cc5] | 1358 | BFS.ShortestPathList[i] = -1;
|
---|
| 1359 | if ((AddedAtomList != NULL) && (AddedAtomList[i] != NULL)) // mark already present atoms (i.e. Root and maybe others) as visited
|
---|
| 1360 | BFS.ColorList[i] = lightgray;
|
---|
| 1361 | else
|
---|
| 1362 | BFS.ColorList[i] = white;
|
---|
| 1363 | }
|
---|
[920c70] | 1364 | //BFS.ShortestPathList[Root->nr] = 0; // done by Calloc
|
---|
[9eefda] | 1365 | }
|
---|
| 1366 | ;
|
---|
[ce7cc5] | 1367 |
|
---|
| 1368 | void BreadthFirstSearchAdd_Free(struct BFSAccounting &BFS)
|
---|
| 1369 | {
|
---|
[920c70] | 1370 | delete[](BFS.PredecessorList);
|
---|
| 1371 | delete[](BFS.ShortestPathList);
|
---|
| 1372 | delete[](BFS.ColorList);
|
---|
[9eefda] | 1373 | delete (BFS.BFSStack);
|
---|
[ce7cc5] | 1374 | BFS.AtomCount = 0;
|
---|
[9eefda] | 1375 | }
|
---|
| 1376 | ;
|
---|
[ce7cc5] | 1377 |
|
---|
[e138de] | 1378 | void BreadthFirstSearchAdd_UnvisitedNode(molecule *Mol, struct BFSAccounting &BFS, atom *&Walker, atom *&OtherAtom, bond *&Binder, bond *&Bond, atom **&AddedAtomList, bond **&AddedBondList, bool IsAngstroem)
|
---|
[ce7cc5] | 1379 | {
|
---|
| 1380 | if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
|
---|
| 1381 | BFS.ColorList[OtherAtom->nr] = lightgray;
|
---|
[9eefda] | 1382 | BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
|
---|
| 1383 | BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
|
---|
[68f03d] | 1384 | DoLog(2) && (Log() << Verbose(2) << "Coloring OtherAtom " << OtherAtom->getName() << " " << ((BFS.ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->getName() << " and its Shortest Path is " << BFS.ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl);
|
---|
[9eefda] | 1385 | if ((((BFS.ShortestPathList[OtherAtom->nr] < BFS.BondOrder) && (Binder != Bond)))) { // Check for maximum distance
|
---|
[a67d19] | 1386 | DoLog(3) && (Log() << Verbose(3));
|
---|
[ce7cc5] | 1387 | if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
|
---|
| 1388 | AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
|
---|
[68f03d] | 1389 | DoLog(0) && (Log() << Verbose(0) << "Added OtherAtom " << OtherAtom->getName());
|
---|
[ce7cc5] | 1390 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
[a67d19] | 1391 | DoLog(0) && (Log() << Verbose(0) << " and bond " << *(AddedBondList[Binder->nr]) << ", ");
|
---|
[9eefda] | 1392 | } else { // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
|
---|
[68f03d] | 1393 | DoLog(0) && (Log() << Verbose(0) << "Not adding OtherAtom " << OtherAtom->getName());
|
---|
[ce7cc5] | 1394 | if (AddedBondList[Binder->nr] == NULL) {
|
---|
| 1395 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
[a67d19] | 1396 | DoLog(0) && (Log() << Verbose(0) << ", added Bond " << *(AddedBondList[Binder->nr]));
|
---|
[ce7cc5] | 1397 | } else
|
---|
[a67d19] | 1398 | DoLog(0) && (Log() << Verbose(0) << ", not added Bond ");
|
---|
[ce7cc5] | 1399 | }
|
---|
[a67d19] | 1400 | DoLog(0) && (Log() << Verbose(0) << ", putting OtherAtom into queue." << endl);
|
---|
[a564be] | 1401 | BFS.BFSStack->push_front(OtherAtom);
|
---|
[ce7cc5] | 1402 | } else { // out of bond order, then replace
|
---|
| 1403 | if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
|
---|
| 1404 | BFS.ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
|
---|
| 1405 | if (Binder == Bond)
|
---|
[a67d19] | 1406 | DoLog(3) && (Log() << Verbose(3) << "Not Queueing, is the Root bond");
|
---|
[ce7cc5] | 1407 | else if (BFS.ShortestPathList[OtherAtom->nr] >= BFS.BondOrder)
|
---|
[a67d19] | 1408 | DoLog(3) && (Log() << Verbose(3) << "Not Queueing, is out of Bond Count of " << BFS.BondOrder);
|
---|
[ce7cc5] | 1409 | if (!Binder->Cyclic)
|
---|
[a67d19] | 1410 | DoLog(0) && (Log() << Verbose(0) << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl);
|
---|
[ce7cc5] | 1411 | if (AddedBondList[Binder->nr] == NULL) {
|
---|
| 1412 | if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
|
---|
| 1413 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
| 1414 | } else {
|
---|
[9eefda] | 1415 | #ifdef ADDHYDROGEN
|
---|
[e138de] | 1416 | if (!Mol->AddHydrogenReplacementAtom(Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
|
---|
[9eefda] | 1417 | exit(1);
|
---|
| 1418 | #endif
|
---|
[ce7cc5] | 1419 | }
|
---|
| 1420 | }
|
---|
| 1421 | }
|
---|
[9eefda] | 1422 | }
|
---|
| 1423 | ;
|
---|
[ce7cc5] | 1424 |
|
---|
[e138de] | 1425 | void BreadthFirstSearchAdd_VisitedNode(molecule *Mol, struct BFSAccounting &BFS, atom *&Walker, atom *&OtherAtom, bond *&Binder, bond *&Bond, atom **&AddedAtomList, bond **&AddedBondList, bool IsAngstroem)
|
---|
[ce7cc5] | 1426 | {
|
---|
[a67d19] | 1427 | DoLog(3) && (Log() << Verbose(3) << "Not Adding, has already been visited." << endl);
|
---|
[ce7cc5] | 1428 | // This has to be a cyclic bond, check whether it's present ...
|
---|
| 1429 | if (AddedBondList[Binder->nr] == NULL) {
|
---|
[9eefda] | 1430 | if ((Binder != Bond) && (Binder->Cyclic) && (((BFS.ShortestPathList[Walker->nr] + 1) < BFS.BondOrder))) {
|
---|
[ce7cc5] | 1431 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
| 1432 | } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
|
---|
[9eefda] | 1433 | #ifdef ADDHYDROGEN
|
---|
[e138de] | 1434 | if(!Mol->AddHydrogenReplacementAtom(Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
|
---|
[9eefda] | 1435 | exit(1);
|
---|
| 1436 | #endif
|
---|
[ce7cc5] | 1437 | }
|
---|
| 1438 | }
|
---|
[9eefda] | 1439 | }
|
---|
| 1440 | ;
|
---|
[cee0b57] | 1441 |
|
---|
| 1442 | /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
|
---|
[a564be] | 1443 | * Gray vertices are always enqueued in an std::deque<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
|
---|
[cee0b57] | 1444 | * white and putting into queue.
|
---|
| 1445 | * \param *out output stream for debugging
|
---|
| 1446 | * \param *Mol Molecule class to add atoms to
|
---|
| 1447 | * \param **AddedAtomList list with added atom pointers, index is atom father's number
|
---|
| 1448 | * \param **AddedBondList list with added bond pointers, index is bond father's number
|
---|
| 1449 | * \param *Root root vertex for BFS
|
---|
| 1450 | * \param *Bond bond not to look beyond
|
---|
| 1451 | * \param BondOrder maximum distance for vertices to add
|
---|
| 1452 | * \param IsAngstroem lengths are in angstroem or bohrradii
|
---|
| 1453 | */
|
---|
[e138de] | 1454 | void molecule::BreadthFirstSearchAdd(molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
|
---|
[cee0b57] | 1455 | {
|
---|
[ce7cc5] | 1456 | struct BFSAccounting BFS;
|
---|
[cee0b57] | 1457 | atom *Walker = NULL, *OtherAtom = NULL;
|
---|
[ce7cc5] | 1458 | bond *Binder = NULL;
|
---|
[cee0b57] | 1459 |
|
---|
| 1460 | // add Root if not done yet
|
---|
[9eefda] | 1461 | if (AddedAtomList[Root->nr] == NULL) // add Root if not yet present
|
---|
[cee0b57] | 1462 | AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
|
---|
| 1463 |
|
---|
[ea7176] | 1464 | BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, getAtomCount(), AddedAtomList);
|
---|
[cee0b57] | 1465 |
|
---|
| 1466 | // and go on ... Queue always contains all lightgray vertices
|
---|
[a564be] | 1467 | while (!BFS.BFSStack->empty()) {
|
---|
[cee0b57] | 1468 | // we have to pop the oldest atom from stack. This keeps the atoms on the stack always of the same ShortestPath distance.
|
---|
| 1469 | // e.g. if current atom is 2, push to end of stack are of length 3, but first all of length 2 would be popped. They again
|
---|
| 1470 | // append length of 3 (their neighbours). Thus on stack we have always atoms of a certain length n at bottom of stack and
|
---|
| 1471 | // followed by n+1 till top of stack.
|
---|
[a564be] | 1472 | Walker = BFS.BFSStack->front(); // pop oldest added
|
---|
| 1473 | BFS.BFSStack->pop_front();
|
---|
[9d83b6] | 1474 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 1475 | DoLog(1) && (Log() << Verbose(1) << "Current Walker is: " << Walker->getName() << ", and has " << ListOfBonds.size() << " bonds." << endl);
|
---|
| 1476 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 1477 | Runner != ListOfBonds.end();
|
---|
| 1478 | ++Runner) {
|
---|
[266237] | 1479 | if ((*Runner) != NULL) { // don't look at bond equal NULL
|
---|
[ce7cc5] | 1480 | Binder = (*Runner);
|
---|
[266237] | 1481 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
[68f03d] | 1482 | DoLog(2) && (Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "." << endl);
|
---|
[ce7cc5] | 1483 | if (BFS.ColorList[OtherAtom->nr] == white) {
|
---|
[e138de] | 1484 | BreadthFirstSearchAdd_UnvisitedNode(Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
|
---|
[cee0b57] | 1485 | } else {
|
---|
[e138de] | 1486 | BreadthFirstSearchAdd_VisitedNode(Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
|
---|
[cee0b57] | 1487 | }
|
---|
| 1488 | }
|
---|
| 1489 | }
|
---|
[ce7cc5] | 1490 | BFS.ColorList[Walker->nr] = black;
|
---|
[68f03d] | 1491 | DoLog(1) && (Log() << Verbose(1) << "Coloring Walker " << Walker->getName() << " black." << endl);
|
---|
[cee0b57] | 1492 | }
|
---|
[ce7cc5] | 1493 | BreadthFirstSearchAdd_Free(BFS);
|
---|
[9eefda] | 1494 | }
|
---|
| 1495 | ;
|
---|
[cee0b57] | 1496 |
|
---|
[266237] | 1497 | /** Adds a bond as a copy to a given one
|
---|
| 1498 | * \param *left leftatom of new bond
|
---|
| 1499 | * \param *right rightatom of new bond
|
---|
| 1500 | * \param *CopyBond rest of fields in bond are copied from this
|
---|
| 1501 | * \return pointer to new bond
|
---|
| 1502 | */
|
---|
| 1503 | bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
|
---|
| 1504 | {
|
---|
| 1505 | bond *Binder = AddBond(left, right, CopyBond->BondDegree);
|
---|
| 1506 | Binder->Cyclic = CopyBond->Cyclic;
|
---|
| 1507 | Binder->Type = CopyBond->Type;
|
---|
| 1508 | return Binder;
|
---|
[9eefda] | 1509 | }
|
---|
| 1510 | ;
|
---|
[266237] | 1511 |
|
---|
[e138de] | 1512 | void BuildInducedSubgraph_Init(atom **&ParentList, int AtomCount)
|
---|
[cee0b57] | 1513 | {
|
---|
| 1514 | // reset parent list
|
---|
[920c70] | 1515 | ParentList = new atom*[AtomCount];
|
---|
| 1516 | for (int i=0;i<AtomCount;i++)
|
---|
| 1517 | ParentList[i] = NULL;
|
---|
[a67d19] | 1518 | DoLog(3) && (Log() << Verbose(3) << "Resetting ParentList." << endl);
|
---|
[9eefda] | 1519 | }
|
---|
| 1520 | ;
|
---|
[cee0b57] | 1521 |
|
---|
[e138de] | 1522 | void BuildInducedSubgraph_FillParentList(const molecule *mol, const molecule *Father, atom **&ParentList)
|
---|
[43587e] | 1523 | {
|
---|
[cee0b57] | 1524 | // fill parent list with sons
|
---|
[a67d19] | 1525 | DoLog(3) && (Log() << Verbose(3) << "Filling Parent List." << endl);
|
---|
[9879f6] | 1526 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
| 1527 | ParentList[(*iter)->father->nr] = (*iter);
|
---|
[cee0b57] | 1528 | // Outputting List for debugging
|
---|
[a7b761b] | 1529 | DoLog(4) && (Log() << Verbose(4) << "Son[" << (*iter)->father->nr << "] of " << (*iter)->father << " is " << ParentList[(*iter)->father->nr] << "." << endl);
|
---|
[cee0b57] | 1530 | }
|
---|
[a7b761b] | 1531 | };
|
---|
[43587e] | 1532 |
|
---|
[e138de] | 1533 | void BuildInducedSubgraph_Finalize(atom **&ParentList)
|
---|
[43587e] | 1534 | {
|
---|
[920c70] | 1535 | delete[](ParentList);
|
---|
[9eefda] | 1536 | }
|
---|
| 1537 | ;
|
---|
[43587e] | 1538 |
|
---|
[e138de] | 1539 | bool BuildInducedSubgraph_CreateBondsFromParent(molecule *mol, const molecule *Father, atom **&ParentList)
|
---|
[43587e] | 1540 | {
|
---|
| 1541 | bool status = true;
|
---|
| 1542 | atom *OtherAtom = NULL;
|
---|
[cee0b57] | 1543 | // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
|
---|
[a67d19] | 1544 | DoLog(3) && (Log() << Verbose(3) << "Creating bonds." << endl);
|
---|
[9879f6] | 1545 | for (molecule::const_iterator iter = Father->begin(); iter != Father->end(); ++iter) {
|
---|
| 1546 | if (ParentList[(*iter)->nr] != NULL) {
|
---|
| 1547 | if (ParentList[(*iter)->nr]->father != (*iter)) {
|
---|
[cee0b57] | 1548 | status = false;
|
---|
| 1549 | } else {
|
---|
[9d83b6] | 1550 | const BondList& ListOfBonds = (*iter)->getListOfBonds();
|
---|
| 1551 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 1552 | Runner != ListOfBonds.end();
|
---|
| 1553 | ++Runner) {
|
---|
[9879f6] | 1554 | OtherAtom = (*Runner)->GetOtherAtom((*iter));
|
---|
[cee0b57] | 1555 | if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
|
---|
[a7b761b] | 1556 | DoLog(4) && (Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[(*iter)->nr]->getName() << " and " << ParentList[OtherAtom->nr]->getName() << "." << endl);
|
---|
[9879f6] | 1557 | mol->AddBond(ParentList[(*iter)->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
|
---|
[cee0b57] | 1558 | }
|
---|
| 1559 | }
|
---|
| 1560 | }
|
---|
| 1561 | }
|
---|
| 1562 | }
|
---|
[43587e] | 1563 | return status;
|
---|
[9eefda] | 1564 | }
|
---|
| 1565 | ;
|
---|
[cee0b57] | 1566 |
|
---|
[43587e] | 1567 | /** Adds bond structure to this molecule from \a Father molecule.
|
---|
| 1568 | * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
|
---|
| 1569 | * with end points present in this molecule, bond is created in this molecule.
|
---|
| 1570 | * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
|
---|
| 1571 | * \param *out output stream for debugging
|
---|
| 1572 | * \param *Father father molecule
|
---|
| 1573 | * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
|
---|
| 1574 | * \todo not checked, not fully working probably
|
---|
| 1575 | */
|
---|
[e138de] | 1576 | bool molecule::BuildInducedSubgraph(const molecule *Father)
|
---|
[43587e] | 1577 | {
|
---|
| 1578 | bool status = true;
|
---|
| 1579 | atom **ParentList = NULL;
|
---|
[a67d19] | 1580 | DoLog(2) && (Log() << Verbose(2) << "Begin of BuildInducedSubgraph." << endl);
|
---|
[ea7176] | 1581 | BuildInducedSubgraph_Init(ParentList, Father->getAtomCount());
|
---|
[e138de] | 1582 | BuildInducedSubgraph_FillParentList(this, Father, ParentList);
|
---|
| 1583 | status = BuildInducedSubgraph_CreateBondsFromParent(this, Father, ParentList);
|
---|
| 1584 | BuildInducedSubgraph_Finalize(ParentList);
|
---|
[a67d19] | 1585 | DoLog(2) && (Log() << Verbose(2) << "End of BuildInducedSubgraph." << endl);
|
---|
[cee0b57] | 1586 | return status;
|
---|
[9eefda] | 1587 | }
|
---|
| 1588 | ;
|
---|
[cee0b57] | 1589 |
|
---|
| 1590 | /** For a given keyset \a *Fragment, checks whether it is connected in the current molecule.
|
---|
| 1591 | * \param *Fragment Keyset of fragment's vertices
|
---|
| 1592 | * \return true - connected, false - disconnected
|
---|
| 1593 | * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
|
---|
| 1594 | */
|
---|
[e138de] | 1595 | bool molecule::CheckForConnectedSubgraph(KeySet *Fragment)
|
---|
[cee0b57] | 1596 | {
|
---|
| 1597 | atom *Walker = NULL, *Walker2 = NULL;
|
---|
| 1598 | bool BondStatus = false;
|
---|
| 1599 | int size;
|
---|
| 1600 |
|
---|
[a67d19] | 1601 | DoLog(1) && (Log() << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl);
|
---|
| 1602 | DoLog(2) && (Log() << Verbose(2) << "Disconnected atom: ");
|
---|
[cee0b57] | 1603 |
|
---|
| 1604 | // count number of atoms in graph
|
---|
| 1605 | size = 0;
|
---|
[9eefda] | 1606 | for (KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
|
---|
[cee0b57] | 1607 | size++;
|
---|
| 1608 | if (size > 1)
|
---|
[9eefda] | 1609 | for (KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++) {
|
---|
[cee0b57] | 1610 | Walker = FindAtom(*runner);
|
---|
| 1611 | BondStatus = false;
|
---|
[9eefda] | 1612 | for (KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
|
---|
[cee0b57] | 1613 | Walker2 = FindAtom(*runners);
|
---|
[9d83b6] | 1614 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
| 1615 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
| 1616 | Runner != ListOfBonds.end();
|
---|
| 1617 | ++Runner) {
|
---|
[266237] | 1618 | if ((*Runner)->GetOtherAtom(Walker) == Walker2) {
|
---|
[cee0b57] | 1619 | BondStatus = true;
|
---|
| 1620 | break;
|
---|
| 1621 | }
|
---|
| 1622 | if (BondStatus)
|
---|
| 1623 | break;
|
---|
| 1624 | }
|
---|
| 1625 | }
|
---|
| 1626 | if (!BondStatus) {
|
---|
[a67d19] | 1627 | DoLog(0) && (Log() << Verbose(0) << (*Walker) << endl);
|
---|
[cee0b57] | 1628 | return false;
|
---|
| 1629 | }
|
---|
| 1630 | }
|
---|
| 1631 | else {
|
---|
[a67d19] | 1632 | DoLog(0) && (Log() << Verbose(0) << "none." << endl);
|
---|
[cee0b57] | 1633 | return true;
|
---|
| 1634 | }
|
---|
[a67d19] | 1635 | DoLog(0) && (Log() << Verbose(0) << "none." << endl);
|
---|
[cee0b57] | 1636 |
|
---|
[a67d19] | 1637 | DoLog(1) && (Log() << Verbose(1) << "End of CheckForConnectedSubgraph" << endl);
|
---|
[cee0b57] | 1638 |
|
---|
| 1639 | return true;
|
---|
| 1640 | }
|
---|