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