source: src/Graph/DepthFirstSearchAnalysis.cpp@ 42127c

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 42127c was 42127c, checked in by Frederik Heber <heber@…>, 13 years ago

Extracted definition of MoleculeListClass and put into own header module.

  • Property mode set to 100644
File size: 16.4 KB
Line 
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
8/*
9 * DepthFirstSearchAnalysis.cpp
10 *
11 * Created on: Feb 16, 2011
12 * Author: heber
13 */
14
15// include config.h
16#ifdef HAVE_CONFIG_H
17#include <config.h>
18#endif
19
20#include "CodePatterns/MemDebug.hpp"
21
22#include "DepthFirstSearchAnalysis.hpp"
23
24#include <algorithm>
25#include <functional>
26
27#include "atom.hpp"
28#include "Bond/bond.hpp"
29#include "CodePatterns/Assert.hpp"
30#include "CodePatterns/Info.hpp"
31#include "CodePatterns/Log.hpp"
32#include "CodePatterns/Verbose.hpp"
33#include "Descriptors/AtomDescriptor.hpp"
34#include "molecule.hpp"
35#include "MoleculeLeafClass.hpp"
36#include "MoleculeListClass.hpp"
37#include "World.hpp"
38
39DepthFirstSearchAnalysis::DepthFirstSearchAnalysis() :
40 CurrentGraphNr(0),
41 ComponentNumber(0),
42 BackStepping(false)
43{
44 ResetAllBondsToUnused();
45}
46
47DepthFirstSearchAnalysis::~DepthFirstSearchAnalysis()
48{}
49
50void DepthFirstSearchAnalysis::Init()
51{
52 CurrentGraphNr = 0;
53 ComponentNumber = 0;
54 BackStepping = false;
55 std::for_each(World::getInstance().getAtomIter(),World::getInstance().atomEnd(),
56 std::mem_fun(&atom::resetGraphNr));
57 std::for_each(World::getInstance().getAtomIter(),World::getInstance().atomEnd(),
58 std::mem_fun(&atom::InitComponentNr));
59}
60
61
62bond * DepthFirstSearchAnalysis::FindNextUnused(atom *vertex) const
63{
64 const BondList& ListOfBonds = vertex->getListOfBonds();
65 for (BondList::const_iterator Runner = ListOfBonds.begin();
66 Runner != ListOfBonds.end();
67 ++Runner)
68 if ((*Runner)->IsUsed() == GraphEdge::white)
69 return ((*Runner));
70 return NULL;
71}
72
73
74void DepthFirstSearchAnalysis::ResetAllBondsToUnused() const
75{
76 World::AtomComposite allatoms = World::getInstance().getAllAtoms();
77 for(World::AtomComposite::const_iterator AtomRunner = allatoms.begin();
78 AtomRunner != allatoms.end();
79 ++AtomRunner) {
80 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
81 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
82 BondRunner != ListOfBonds.end();
83 ++BondRunner)
84 if ((*BondRunner)->leftatom == *AtomRunner)
85 (*BondRunner)->ResetUsed();
86 }
87}
88
89void DepthFirstSearchAnalysis::SetNextComponentNumber(atom *vertex, int nr) const
90{
91 size_t i = 0;
92 ASSERT(vertex != NULL,
93 "DepthFirstSearchAnalysis::SetNextComponentNumber() - Given vertex is NULL!");
94 const BondList& ListOfBonds = vertex->getListOfBonds();
95 for (; i < ListOfBonds.size(); i++) {
96 if (vertex->ComponentNr[i] == -1) { // check if not yet used
97 vertex->ComponentNr[i] = nr;
98 break;
99 } else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
100 break; // breaking here will not cause error!
101 }
102 ASSERT(i < ListOfBonds.size(),
103 "DepthFirstSearchAnalysis::SetNextComponentNumber() - All Component entries are already occupied!");
104}
105
106
107bool DepthFirstSearchAnalysis::PickLocalBackEdges(atom **ListOfLocalAtoms, std::deque<bond *> *&LocalStack) const
108{
109 bool status = true;
110 if (BackEdgeStack.empty()) {
111 ELOG(1, "Reference BackEdgeStack is empty!");
112 return false;
113 }
114 bond *Binder = BackEdgeStack.front();
115 bond *FirstBond = Binder; // mark the first bond, so that we don't loop through the stack indefinitely
116 atom *Walker = NULL, *OtherAtom = NULL;
117
118 do { // go through all bonds and push local ones
119 Walker = ListOfLocalAtoms[Binder->leftatom->getNr()]; // get one atom in the reference molecule
120 if (Walker != NULL) { // if this Walker exists in the subgraph ...
121 const BondList& ListOfBonds = Walker->getListOfBonds();
122 for (BondList::const_iterator Runner = ListOfBonds.begin();
123 Runner != ListOfBonds.end();
124 ++Runner) {
125 OtherAtom = (*Runner)->GetOtherAtom(Walker);
126 if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->getNr()]) { // found the bond
127 LocalStack->push_front((*Runner));
128 LOG(3, "INFO: Found local edge " << *(*Runner) << ".");
129 break;
130 }
131 }
132 }
133 ASSERT(!BackEdgeStack.empty(), "DepthFirstSearchAnalysis::PickLocalBackEdges() - ReferenceStack is empty!");
134 Binder = BackEdgeStack.front(); // loop the stack for next item
135 LOG(3, "Current candidate edge " << Binder << ".");
136 } while (FirstBond != Binder);
137
138 return status;
139}
140
141
142
143void DepthFirstSearchAnalysis::OutputGraphInfoPerAtom() const
144{
145 LOG(1, "Final graph info for each atom is:");
146 World::AtomComposite allatoms = World::getInstance().getAllAtoms();
147 for_each(allatoms.begin(),allatoms.end(),mem_fun(&atom::OutputGraphInfo));
148}
149
150
151void DepthFirstSearchAnalysis::OutputGraphInfoPerBond() const
152{
153 LOG(1, "Final graph info for each bond is:");
154 World::AtomComposite allatoms = World::getInstance().getAllAtoms();
155 for(World::AtomComposite::const_iterator AtomRunner = allatoms.begin();
156 AtomRunner != allatoms.end();
157 ++AtomRunner) {
158 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
159 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
160 BondRunner != ListOfBonds.end();
161 ++BondRunner)
162 if ((*BondRunner)->leftatom == *AtomRunner) {
163 const bond *Binder = *BondRunner;
164 if (DoLog(2)) {
165 ostream &out = (Log() << Verbose(2));
166 out << ((Binder->Type == GraphEdge::TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
167 out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
168 Binder->leftatom->OutputComponentNumber(&out);
169 out << " === ";
170 out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
171 Binder->rightatom->OutputComponentNumber(&out);
172 out << ">." << endl;
173 }
174 if (Binder->Cyclic) // cyclic ??
175 LOG(3, "Lowpoint at each side are equal: CYCLIC!");
176 }
177 }
178}
179
180
181unsigned int DepthFirstSearchAnalysis::CyclicBondAnalysis() const
182{
183 unsigned int NoCyclicBonds = 0;
184 World::AtomComposite allatoms = World::getInstance().getAllAtoms();
185 for(World::AtomComposite::const_iterator AtomRunner = allatoms.begin();
186 AtomRunner != allatoms.end();
187 ++AtomRunner) {
188 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
189 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
190 BondRunner != ListOfBonds.end();
191 ++BondRunner)
192 if ((*BondRunner)->leftatom == *AtomRunner)
193 if ((*BondRunner)->rightatom->LowpointNr == (*BondRunner)->leftatom->LowpointNr) { // cyclic ??
194 (*BondRunner)->Cyclic = true;
195 NoCyclicBonds++;
196 }
197 }
198 return NoCyclicBonds;
199}
200
201
202void DepthFirstSearchAnalysis::SetWalkersGraphNr(atom *&Walker)
203{
204 if (!BackStepping) { // if we don't just return from (8)
205 Walker->GraphNr = CurrentGraphNr;
206 Walker->LowpointNr = CurrentGraphNr;
207 LOG(1, "Setting Walker[" << Walker->getName() << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << ".");
208 AtomStack.push_front(Walker);
209 CurrentGraphNr++;
210 }
211}
212
213
214void DepthFirstSearchAnalysis::ProbeAlongUnusedBond(atom *&Walker, bond *&Binder)
215{
216 atom *OtherAtom = NULL;
217
218 do { // (3) if Walker has no unused egdes, go to (5)
219 BackStepping = false; // reset backstepping flag for (8)
220 if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
221 Binder = FindNextUnused(Walker);
222 if (Binder == NULL)
223 break;
224 LOG(2, "Current Unused Bond is " << *Binder << ".");
225 // (4) Mark Binder used, ...
226 Binder->MarkUsed(GraphEdge::black);
227 OtherAtom = Binder->GetOtherAtom(Walker);
228 LOG(2, "(4) OtherAtom is " << OtherAtom->getName() << ".");
229 if (OtherAtom->GraphNr != -1) {
230 // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
231 Binder->Type = GraphEdge::BackEdge;
232 BackEdgeStack.push_front(Binder);
233 Walker->LowpointNr = (Walker->LowpointNr < OtherAtom->GraphNr) ? Walker->LowpointNr : OtherAtom->GraphNr;
234 LOG(3, "(4a) Visited: Setting Lowpoint of Walker[" << Walker->getName() << "] to " << Walker->LowpointNr << ".");
235 } else {
236 // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
237 Binder->Type = GraphEdge::TreeEdge;
238 OtherAtom->Ancestor = Walker;
239 Walker = OtherAtom;
240 LOG(3, "(4b) Not Visited: OtherAtom[" << OtherAtom->getName() << "]'s Ancestor is now " << OtherAtom->Ancestor->getName() << ", Walker is OtherAtom " << OtherAtom->getName() << ".");
241 break;
242 }
243 Binder = NULL;
244 } while (1); // (3)
245}
246
247
248void DepthFirstSearchAnalysis::CheckForaNewComponent(atom *&Walker, ConnectedSubgraph &Subgraph)
249{
250 atom *OtherAtom = NULL;
251
252 // (5) if Ancestor of Walker is ...
253 LOG(1, "(5) Number of Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "] is " << Walker->Ancestor->GraphNr << ".");
254
255 if (Walker->Ancestor->GraphNr != Root->GraphNr) {
256 // (6) (Ancestor of Walker is not Root)
257 if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
258 // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
259 Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
260 LOG(2, "(6) Setting Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << ".");
261 } else {
262 // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
263 Walker->Ancestor->SeparationVertex = true;
264 LOG(2, "(7) Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s is a separating vertex, creating component.");
265 SetNextComponentNumber(Walker->Ancestor, ComponentNumber);
266 LOG(3, "(7) Walker[" << Walker->getName() << "]'s Ancestor's Compont is " << ComponentNumber << ".");
267 SetNextComponentNumber(Walker, ComponentNumber);
268 LOG(3, "(7) Walker[" << Walker->getName() << "]'s Compont is " << ComponentNumber << ".");
269 do {
270 ASSERT(!AtomStack.empty(), "DepthFirstSearchAnalysis_CheckForaNewComponent() - AtomStack is empty!");
271 OtherAtom = AtomStack.front();
272 AtomStack.pop_front();
273 Subgraph.push_back(OtherAtom);
274 SetNextComponentNumber(OtherAtom, ComponentNumber);
275 LOG(3, "(7) Other[" << OtherAtom->getName() << "]'s Compont is " << ComponentNumber << ".");
276 } while (OtherAtom != Walker);
277 ComponentNumber++;
278 }
279 // (8) Walker becomes its Ancestor, go to (3)
280 LOG(2, "(8) Walker[" << Walker->getName() << "] is now its Ancestor " << Walker->Ancestor->getName() << ", backstepping. ");
281 Walker = Walker->Ancestor;
282 BackStepping = true;
283 }
284}
285
286
287void DepthFirstSearchAnalysis::CleanRootStackDownTillWalker(atom *&Walker, bond *&Binder, ConnectedSubgraph &Subgraph)
288{
289 atom *OtherAtom = NULL;
290
291 if (!BackStepping) { // coming from (8) want to go to (3)
292 // (9) remove all from stack till Walker (including), these and Root form a component
293 //AtomStack.Output(out);
294 SetNextComponentNumber(Root, ComponentNumber);
295 LOG(3, "(9) Root[" << Root->getName() << "]'s Component is " << ComponentNumber << ".");
296 SetNextComponentNumber(Walker, ComponentNumber);
297 LOG(3, "(9) Walker[" << Walker->getName() << "]'s Component is " << ComponentNumber << ".");
298 do {
299 ASSERT(!AtomStack.empty(), "DepthFirstSearchAnalysis::CleanRootStackDownTillWalker() - AtomStack is empty!");
300 OtherAtom = AtomStack.front();
301 AtomStack.pop_front();
302 Subgraph.push_back(OtherAtom);
303 SetNextComponentNumber(OtherAtom, ComponentNumber);
304 LOG(3, "(7) Other[" << OtherAtom->getName() << "]'s Component is " << ComponentNumber << ".");
305 } while (OtherAtom != Walker);
306 ComponentNumber++;
307
308 // (11) Root is separation vertex, set Walker to Root and go to (4)
309 Walker = Root;
310 Binder = FindNextUnused(Walker);
311 if (Binder != NULL) { // Root is separation vertex
312 LOG(1, "(10) Walker is Root[" << Root->getName() << "], next Unused Bond is " << *Binder << ".");
313 LOG(1, "(11) Root is a separation vertex.");
314 Walker->SeparationVertex = true;
315 } else {
316 LOG(1, "(10) Walker is Root[" << Root->getName() << "], no next Unused Bond.");
317 }
318 }
319}
320
321
322const std::deque<bond *>& DepthFirstSearchAnalysis::getBackEdgeStack() const
323{
324 return BackEdgeStack;
325}
326
327
328void DepthFirstSearchAnalysis::operator()()
329{
330 Info FunctionInfo("DepthFirstSearchAnalysis");
331 ListOfConnectedSubgraphs.clear();
332 int OldGraphNr = 0;
333 atom *Walker = NULL;
334 bond *Binder = NULL;
335
336 if (World::getInstance().numAtoms() == 0)
337 return;
338
339 Init();
340
341 LOG(0, "STATUS: Start walking the bond graph.");
342 for(World::AtomIterator iter = World::getInstance().getAtomIter();
343 iter != World::getInstance().atomEnd();) { // don't advance, is done at the end
344 Root = *iter;
345 // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all
346 AtomStack.clear();
347
348 // put into new subgraph molecule and add this to list of subgraphs
349 ConnectedSubgraph CurrentSubgraph;
350 CurrentSubgraph.push_back(Root);
351
352 OldGraphNr = CurrentGraphNr;
353 Walker = Root;
354 do { // (10)
355 do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom
356 SetWalkersGraphNr(Walker);
357
358 ProbeAlongUnusedBond(Walker, Binder);
359
360 if (Binder == NULL) {
361 LOG(2, "No more Unused Bonds.");
362 break;
363 } else
364 Binder = NULL;
365 } while (1); // (2)
366
367 // 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!
368 if ((Walker == Root) && (Binder == NULL))
369 break;
370
371 CheckForaNewComponent( Walker, CurrentSubgraph);
372
373 CleanRootStackDownTillWalker(Walker, Binder, CurrentSubgraph);
374
375 } while ((BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
376
377 ListOfConnectedSubgraphs.push_back(CurrentSubgraph);
378 // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
379 std::stringstream output;
380 output << CurrentSubgraph;
381 LOG(0, "STATUS: Disconnected subgraph ranges from " << OldGraphNr << " to "
382 << CurrentGraphNr-1 << ": " << output.str());
383
384 // step on to next root
385 while (iter != World::getInstance().atomEnd()) {
386 if ((*iter)->GraphNr != -1) { // if already discovered, step on
387 iter++;
388 } else {
389 LOG(1,"Current next subgraph root candidate is " << (*iter)->getName()
390 << " with GraphNr " << (*iter)->GraphNr << ".");
391 break;
392 }
393 }
394 }
395 LOG(0, "STATUS: Done walking the bond graph.");
396
397 // set cyclic bond criterium on "same LP" basis
398 CyclicBondAnalysis();
399
400 OutputGraphInfoPerAtom();
401
402 OutputGraphInfoPerBond();
403}
404
405void DepthFirstSearchAnalysis::UpdateMoleculeStructure() const
406{
407 // remove all of World's molecules
408 for (World::MoleculeIterator iter = World::getInstance().getMoleculeIter();
409 World::getInstance().getMoleculeIter() != World::getInstance().moleculeEnd();
410 iter = World::getInstance().getMoleculeIter()) {
411 World::getInstance().getMolecules()->erase(*iter);
412 World::getInstance().destroyMolecule(*iter);
413 }
414 // instantiate new molecules
415 molecule *newmol = NULL;
416 for (ConnectedSubgraphList::const_iterator iter = ListOfConnectedSubgraphs.begin();
417 iter != ListOfConnectedSubgraphs.end();
418 ++iter) {
419 LOG(0, "STATUS: Creating new molecule:");
420 std::stringstream output;
421 newmol = (*iter).getMolecule();
422 newmol->Output(&output);
423 std::stringstream outstream(output.str());
424 std::string line;
425 while (getline(outstream, line)) {
426 LOG(0, "\t"+line);
427 }
428 }
429}
430
431MoleculeLeafClass *DepthFirstSearchAnalysis::getMoleculeStructure() const
432{
433 MoleculeLeafClass *Subgraphs = new MoleculeLeafClass(NULL);
434 MoleculeLeafClass *MolecularWalker = Subgraphs;
435 for (World::MoleculeIterator iter = World::getInstance().getMoleculeIter();
436 iter != World::getInstance().moleculeEnd();
437 ++iter) {
438 // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
439 MolecularWalker = new MoleculeLeafClass(MolecularWalker);
440 MolecularWalker->Leaf = (*iter);
441 }
442 return Subgraphs;
443}
444
Note: See TracBrowser for help on using the repository browser.