source: src/FunctionApproximation/Extractors.cpp@ c099dc

Candidate_v1.7.0 stable
Last change on this file since c099dc was d951a5, checked in by Frederik Heber <frederik.heber@…>, 5 years ago

FIX: Small changes to ChemicalSpaceEvaluatorAction.

  • discarding degree combination before we begin calculating the fragment

energy.

  • made the log info more specific and easier to read.
  • FIX: CSE used wrong graph when generating the homology graph key.
  • Property mode set to 100644
File size: 32.7 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2012 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
7 *
8 *
9 * This file is part of MoleCuilder.
10 *
11 * MoleCuilder is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * MoleCuilder is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
23 */
24
25/*
26 * Extractors.cpp
27 *
28 * Created on: 15.10.2012
29 * Author: heber
30 */
31
32// include config.h
33#ifdef HAVE_CONFIG_H
34#include <config.h>
35#endif
36
37//#include "CodePatterns/MemDebug.hpp"
38
39#include <algorithm>
40#include <iterator>
41#include <set>
42#include <sstream>
43#include <utility>
44#include <vector>
45#include <boost/assign.hpp>
46#include <boost/bind.hpp>
47#include <boost/foreach.hpp>
48
49#include "CodePatterns/Assert.hpp"
50#include "CodePatterns/IteratorAdaptors.hpp"
51#include "CodePatterns/Log.hpp"
52#include "CodePatterns/toString.hpp"
53
54#include "LinearAlgebra/Vector.hpp"
55
56#include "Fragmentation/Homology/HomologyGraph.hpp"
57#include "FunctionApproximation/Extractors.hpp"
58#include "FunctionApproximation/FunctionArgument.hpp"
59#include "Potentials/BindingModel.hpp"
60
61using namespace boost::assign;
62
63using namespace Extractors;
64
65FunctionModel::arguments_t
66Extractors::gatherAllSymmetricDistanceArguments(
67 const Fragment::positions_t& positions,
68 const Fragment::atomicnumbers_t& atomicnumbers,
69 const FragmentationEdges::edges_t &edges,
70 const size_t globalid)
71{
72 FunctionModel::arguments_t result;
73
74 // place edges in map
75 typedef std::set< std::pair<size_t, size_t> > sorted_edges_t;
76 sorted_edges_t sorted_edges;
77 for (FragmentationEdges::edges_t::const_iterator edgeiter = edges.begin();
78 edgeiter != edges.end(); ++edgeiter) {
79 std::pair<sorted_edges_t::iterator, bool> inserter =
80 sorted_edges.insert(
81 (edgeiter->first < edgeiter->second) ?
82 std::make_pair(edgeiter->first, edgeiter->second) :
83 std::make_pair(edgeiter->second, edgeiter->first));
84 ASSERT(inserter.second,
85 "Extractors::gatherAllSymmetricDistanceArguments() - edge "
86 +toString(*inserter.first)+" is already present");
87 }
88
89 // go through current configuration and gather all other distances
90 Fragment::positions_t::const_iterator firstpositer = positions.begin();
91 for (;firstpositer != positions.end(); ++firstpositer) {
92 Fragment::positions_t::const_iterator secondpositer = firstpositer;
93 for (; secondpositer != positions.end(); ++secondpositer) {
94 if (firstpositer == secondpositer)
95 continue;
96 argument_t arg;
97 const Vector firsttemp((*firstpositer)[0],(*firstpositer)[1],(*firstpositer)[2]);
98 const Vector secondtemp((*secondpositer)[0],(*secondpositer)[1],(*secondpositer)[2]);
99 arg.distance = firsttemp.distance(secondtemp);
100 arg.types = std::make_pair(
101 (int)atomicnumbers[ std::distance(positions.begin(), firstpositer) ],
102 (int)atomicnumbers[ std::distance(positions.begin(), secondpositer) ]
103 );
104 arg.indices = std::make_pair(
105 std::distance(
106 positions.begin(), firstpositer),
107 std::distance(
108 positions.begin(), secondpositer)
109 );
110 arg.globalid = globalid;
111 arg.bonded = sorted_edges.find(arg.indices) != sorted_edges.end();
112 LOG(3, "DEBUG: Created argument " << arg << ".");
113 result.push_back(arg);
114 }
115 }
116
117 return result;
118}
119
120Extractors::elementcounts_t
121Extractors::_detail::getElementCounts(
122 const Fragment::atomicnumbers_t elements
123 )
124{
125 elementcounts_t elementcounts;
126 for (Fragment::atomicnumbers_t::const_iterator elementiter = elements.begin();
127 elementiter != elements.end(); ++elementiter) {
128 // insert new element
129 std::pair< elementcounts_t::iterator, bool> inserter =
130 elementcounts.insert( std::make_pair( *elementiter, 1) );
131 // if already present, just increase its count
132 if (!inserter.second)
133 ++(inserter.first->second);
134 }
135 return elementcounts;
136}
137
138struct ParticleTypesComparator {
139 bool operator()(const argument_t::types_t &a, const argument_t::types_t &b)
140 {
141 if (a.first < a.second) {
142 if (b.first < b.second) {
143 if (a.first < b.first)
144 return true;
145 else if (a.first > b.first)
146 return false;
147 else
148 return (a.second < b.second);
149 } else {
150 if (a.first < b.second)
151 return true;
152 else if (a.first > b.second)
153 return false;
154 else
155 return (a.second < b.first);
156 }
157 } else {
158 if (b.first < b.second) {
159 if (a.second < b.first)
160 return true;
161 else if (a.second > b.first)
162 return false;
163 else
164 return (a.first < b.second);
165 } else {
166 if (a.second < b.second)
167 return true;
168 else if (a.second > b.second)
169 return false;
170 else
171 return (a.first < b.first);
172 }
173 }
174 }
175};
176
177std::ostream& operator<<(std::ostream &out, const argument_t::types_t &a)
178{
179 out << "[" << a.first << "," << a.second << "]";
180 return out;
181}
182
183void insertIntoNodeFragmentMap(
184 node_FragmentNode_map_t &_node_FragmentNode_map,
185 const size_t &_index,
186 const Extractors::ParticleType_t &_type)
187{
188 const node_FragmentNode_map_t::iterator mapiter = _node_FragmentNode_map.find(_index);
189 // check if already present
190 if (mapiter != _node_FragmentNode_map.end()) {
191 // assert same type and increment number of edges
192 ASSERT( mapiter->second.first == _type,
193 "insertIntoNodeFragmentMap() - different types "+toString(mapiter->second.first)+
194 " and "+toString(_type)+" for node "+toString(_index));
195 ++(mapiter->second.second);
196 } else {
197 // place new entry with a single present edge
198 _node_FragmentNode_map.insert( std::make_pair(_index, std::make_pair(_type, 1) ));
199 }
200}
201
202static node_FragmentNode_map_t fillNodeFragmentMap(
203 FunctionModel::arguments_t &argumentbunch)
204{
205 node_FragmentNode_map_t node_FragmentNode_map;
206 // place each node index with type and number of edges into map
207 for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
208 argiter != argumentbunch.end(); ++argiter) {
209 const argument_t &arg = *argiter;
210 // only consider the distances that represent a bond edge
211 if (arg.bonded) {
212 insertIntoNodeFragmentMap(node_FragmentNode_map, arg.indices.first, arg.types.first);
213 insertIntoNodeFragmentMap(node_FragmentNode_map, arg.indices.second, arg.types.second);
214 }
215 }
216
217 return node_FragmentNode_map;
218}
219
220static argument_placement_map_t fillArgumentsPlacementMap(const size_t num_args)
221{
222 argument_placement_map_t argument_placement_map;
223 size_t placement = 0;
224 for (size_t i = 0;i<num_args;++i)
225 for (size_t j = i+1;j<num_args;++j)
226 argument_placement_map.insert( std::make_pair( argument_t::indices_t(i,j), placement++) );
227 return argument_placement_map;
228}
229
230static argument_t::indices_t translateIndices(
231 const argindex_to_nodeindex_t &_argindex_to_nodeindex,
232 const argument_t::indices_t &_indices
233 )
234{
235 const argindex_to_nodeindex_t::const_iterator firstiter =
236 _argindex_to_nodeindex.find(_indices.first);
237 ASSERT( firstiter != _argindex_to_nodeindex.end(),
238 "translateIndices() - could not find #"+toString(_indices.first)+
239 " in translation map.");
240 const argindex_to_nodeindex_t::const_iterator seconditer =
241 _argindex_to_nodeindex.find(_indices.second);
242 ASSERT(seconditer != _argindex_to_nodeindex.end(),
243 "translateIndices() - could not find #"+toString(_indices.second)+
244 " in translation map.");
245 // mind increasing order of indices in pair
246 if (firstiter->second < seconditer->second)
247 return argument_t::indices_t(firstiter->second, seconditer->second);
248 else
249 return argument_t::indices_t(seconditer->second, firstiter->second);
250}
251
252/** Power set generator
253 *
254 * taken from https://rosettacode.org/wiki/Power_set#Non-recursive_version
255 *
256 */
257static powerset_type powerset(set_type const& set)
258{
259 typedef set_type::const_iterator set_iter;
260 typedef std::vector<set_iter> vec;
261
262 struct local
263 {
264 static int dereference(set_iter v) { return *v; }
265 };
266
267 powerset_type result;
268
269 vec elements;
270 do {
271 set_type tmp;
272 std::transform(elements.begin(), elements.end(),
273 std::inserter(tmp, tmp.end()),
274 local::dereference);
275 result.insert(tmp);
276 if (!elements.empty() && ++elements.back() == set.end()) {
277 elements.pop_back();
278 } else {
279 set_iter iter;
280 if (elements.empty()) {
281 iter = set.begin();
282 } else {
283 iter = elements.back();
284 ++iter;
285 }
286 for (; iter != set.end(); ++iter) {
287 elements.push_back(iter);
288 }
289 }
290 } while (!elements.empty());
291
292 return result;
293}
294
295/** Recursive function to generate all induced, connected subgraphs given a
296 * graph.
297 *
298 * \param N number of still left to pick
299 * \param depth level in \a set_of_nodes
300 * \param nodes current set of nodes that are picked already
301 * \param set_of_nodes resulting set of generated subgraphs' nodes
302 * \param nodes_per_level level-wise frontier of connected nodes around a root node
303 * \param graph graph containing the adjacency
304 * \param index_map with indices per \a graph' vertex
305 */
306void Extractors::generateAllInducedConnectedSubgraphs(
307 const size_t N,
308 const level_t level,
309 const nodes_t &nodes,
310 set_of_nodes_t &set_of_nodes,
311 const nodes_per_level_t &nodes_per_level,
312 const UndirectedGraph &graph,
313 const std::vector<size_t> &_distance,
314 const index_map_t &index_map)
315{
316 ASSERT( nodes_per_level.find(level) != nodes_per_level.end(),
317 "generateAllInducedConnectedSubgraphs() - we are deeper than the graph.");
318 ASSERT( N < nodes_per_level.size(),
319 "generateAllInducedConnectedSubgraphs() - we are looking for subgraphs larger than the graph.");
320 if (N > 0) {
321 LOG(3, "DEBUG: At level " << level << " current nodes is " << nodes << ", need to find " << N << " more.");
322 // get next level's set and constrain to nodes connected to this set
323 nodes_t validnodes;
324 std::pair< nodes_per_level_t::const_iterator, nodes_per_level_t::const_iterator> range =
325 nodes_per_level.equal_range(level);
326 for (nodes_per_level_t::const_iterator rangeiter = range.first;
327 rangeiter != range.second; ++rangeiter) {
328 LOG(4, "DEBUG: Checking edges further away from node #" << rangeiter->second);
329 // get all edges connected to this node further away
330 UndirectedGraph::in_edge_iterator i, end;
331 boost::tie(i, end) = boost::in_edges(boost::vertex(rangeiter->second, graph), graph);
332 for (;i != end; ++i) {
333 // check each edge whether it's in nodes
334 const node_t sourceindex = boost::get(index_map, boost::source(*i, graph));
335 const node_t targetindex = boost::get(index_map, boost::target(*i, graph));
336 const size_t &source_distance = _distance[sourceindex];
337 const size_t &target_distance = _distance[targetindex];
338 // edge is going deeper into graph
339 if (((source_distance == level) && (target_distance == (level+1)))
340 || ((source_distance == (level+1)) && (target_distance == level))) {
341 LOG(5, "DEBUG: Candidate edge on level " << level << " is from " << sourceindex
342 << " to " << targetindex << ".");
343 // pick right index and check for not present in list yet
344 if (sourceindex == rangeiter->second) {
345 if (nodes.count(targetindex) == 0) {
346 validnodes.insert(targetindex);
347 LOG(4, "DEBUG: Inserting node #" << targetindex << " into valid nodes.");
348 }
349 } else if (targetindex == rangeiter->second) {
350 if (nodes.count(sourceindex) == 0) {
351 validnodes.insert(sourceindex);
352 LOG(4, "DEBUG: Inserting node #" << sourceindex << " into valid nodes.");
353 }
354 } else {
355 ASSERT(0,
356 "generateAllInducedConnectedSubgraphs() - neither source #"+toString(sourceindex)+
357 " nor target #"+toString(targetindex)+" is equal to #"+toString(rangeiter->second));
358 }
359 }
360 }
361 }
362 // skip this if we cannot go deeper into the graph from here
363 if (validnodes.empty()) {
364 LOG(3, "DEBUG: We did not find any more nodes to step on from " << nodes << ".");
365 return;
366 }
367
368 // go through power set
369 const powerset_type test_powerset = powerset(validnodes);
370 for (powerset_type::const_iterator iter = test_powerset.begin();
371 iter != test_powerset.end();
372 ++iter) {
373 // count set bits (#elements in *iter), must be between 1 and N
374 const size_t num_set_bits = iter->size();
375 if ((num_set_bits > 0) && (num_set_bits <= N)) {
376 // add set to nodes
377 LOG(3, "DEBUG: With present " << nodes << " the current set is " << *iter << " of "
378 << validnodes << ".");
379
380 // copy the nodes before insertion
381 nodes_t filled_nodes(nodes.begin(), nodes.end());
382 filled_nodes.insert(iter->begin(), iter->end());
383
384 // and recurse
385 generateAllInducedConnectedSubgraphs(
386 N-num_set_bits, level+1, filled_nodes, set_of_nodes, nodes_per_level, graph, _distance, index_map);
387 }
388 }
389 } else {
390 // N==0: we have a winner
391 std::pair<set_of_nodes_t::iterator, bool> inserter =
392 set_of_nodes.insert( nodes );
393 if (!inserter.second)
394 LOG(2, "DEBUG: subgraph " << nodes << " is already contained in set_of_nodes.");
395 else
396 LOG(2, "DEBUG: subgraph " << nodes << " inserted into set_of_nodes.");
397 }
398}
399
400Extractors::ParticleType_t Extractors::getParticleTypeToNode(
401 const type_index_lookup_t &type_index_lookup,
402 const size_t nodeindex)
403{
404 const type_index_lookup_t::left_const_iterator typeiter = type_index_lookup.left.find(nodeindex);
405 ASSERT( typeiter != type_index_lookup.left.end(),
406 "Extractors::getParticleTypeToNode() - could not find type to node #"+toString(nodeindex));
407 return typeiter->second;
408}
409
410HomologyGraph Extractors::createHomologyGraphFromNodes(
411 const nodes_t &nodes,
412 const type_index_lookup_t &type_index_lookup,
413 const UndirectedGraph &graph,
414 const index_map_t &index_map
415 )
416{
417 HomologyGraph::nodes_t graph_nodes;
418 HomologyGraph::edges_t graph_edges;
419 {
420 typedef std::set< std::pair<node_t, node_t> > graph_edges_t;
421 graph_edges_t edge_set;
422 std::pair<HomologyGraph::nodes_t::iterator, bool > inserter;
423 for (nodes_t::const_iterator nodeiter = nodes.begin();
424 nodeiter != nodes.end(); ++nodeiter) {
425 const Extractors::ParticleType_t &nodetype = getParticleTypeToNode(type_index_lookup, *nodeiter);
426
427 // count edges in constrained set for this particular node
428 size_t num_edges = 0;
429 UndirectedGraph::in_edge_iterator i, end;
430 for (boost::tie(i, end) = boost::in_edges(boost::vertex(*nodeiter, graph), graph);
431 i != end; ++i) {
432 const node_t sourceindex = boost::get(index_map, boost::source(*i, graph));
433 const node_t targetindex = boost::get(index_map, boost::target(*i, graph));
434 // check each edge whether it's in nodes
435 if ((nodes.count(sourceindex) != 0) && (nodes.count(targetindex) != 0)) {
436 // increase edge count of node
437 ++num_edges;
438 // we first store edges in a set to ensure their uniqueness (as we encounter
439 // each edge two times and we don't know if source and target will be
440 // different the second time encountered).
441 LOG(4, "DEBUG: Adding edge " << sourceindex << "<->" << targetindex << " to set.");
442 if (sourceindex < targetindex)
443 edge_set.insert( std::make_pair(sourceindex, targetindex) );
444 else
445 edge_set.insert( std::make_pair(targetindex, sourceindex) );
446 }
447 }
448 LOG(4, "DEBUG: Node #" << *nodeiter << " has " << num_edges << " edges.");
449
450 // add FragmentNode
451 inserter = graph_nodes.insert( std::make_pair(FragmentNode(nodetype, num_edges), 1) );
452 if (!inserter.second)
453 ++(inserter.first->second);
454 }
455
456 // add edges
457 for (graph_edges_t::const_iterator edgeiter = edge_set.begin();
458 edgeiter != edge_set.end(); ++edgeiter) {
459 const Extractors::ParticleType_t sourcetype =
460 getParticleTypeToNode(type_index_lookup, edgeiter->first);
461 const Extractors::ParticleType_t targettype =
462 getParticleTypeToNode(type_index_lookup, edgeiter->second);
463 // FragmentEdge takes care of proper sorting
464 FragmentEdge edge(sourcetype, targettype);
465 LOG(4, "DEBUG: Adding fragment edge " << edge);
466 std::pair<HomologyGraph::edges_t::iterator, bool > inserter;
467 inserter = graph_edges.insert( std::make_pair( edge, 1) );
468 if (!inserter.second)
469 ++(inserter.first->second);
470 }
471 }
472
473 return HomologyGraph(graph_nodes, graph_edges);
474}
475
476FunctionModel::list_of_arguments_t Extractors::filterArgumentsByBindingModel(
477 const FunctionModel::arguments_t &args,
478 const HomologyGraph &_graph,
479 const ParticleTypes_t &_types,
480 const BindingModel &_bindingmodel
481 )
482{
483 FunctionModel::list_of_arguments_t returnargs;
484
485 // deal with the case when there are no distances (ConstantPotential)
486 if (_bindingmodel.getNodes().size() < 2) {
487 LOG(3, "DEBUG: Potential requires only one or no particle types, needs no distances.");
488 return returnargs;
489 }
490 if (_bindingmodel.getGraph().getEdges().empty()) {
491 LOG(3, "DEBUG: Potential represents non-bonded interactions, gets all distances.");
492 // TODO: Here we need to constrain to all distances matching the types?
493 returnargs.push_back(args);
494 return returnargs;
495 }
496
497 /// 0. place all nodes in a lookup map for their type
498 type_index_lookup_t type_index_lookup;
499 for(FunctionModel::arguments_t::const_iterator iter = args.begin();
500 iter != args.end(); ++iter) {
501 if (type_index_lookup.left.find(iter->indices.first) == type_index_lookup.left.end())
502 type_index_lookup.left.insert( std::make_pair(iter->indices.first, iter->types.first) );
503 else {
504 ASSERT(type_index_lookup.left.at(iter->indices.first) == iter->types.first,
505 "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.first)+
506 " is already present with different type "+toString(iter->types.first)
507 +" than "+toString(type_index_lookup.left.at(iter->indices.first)));
508 }
509 if (type_index_lookup.left.find(iter->indices.second) == type_index_lookup.left.end())
510 type_index_lookup.left.insert( std::make_pair(iter->indices.second, iter->types.second) );
511 else {
512 ASSERT(type_index_lookup.left.at(iter->indices.second) == iter->types.second,
513 "Extractors::reorderArgumentsByParticleTypes() - entry " +toString(iter->indices.second)+
514 " is already present with different type "+toString(iter->types.second)
515 +" than "+toString(type_index_lookup.left.at(iter->indices.second)));
516 }
517 }
518 if (DoLog(3)) {
519 std::stringstream output;
520 for (type_index_lookup_t::left_const_iterator indexiter = type_index_lookup.left.begin();
521 indexiter != type_index_lookup.left.end(); ++indexiter) {
522 output << " [" << indexiter->first << "," << indexiter->second << "]";
523 }
524 LOG(3, "DEBUG: index to type map:" << output.str());
525 }
526 if (DoLog(3)) {
527 std::stringstream output;
528 for (type_index_lookup_t::right_const_iterator typeiter = type_index_lookup.right.begin();
529 typeiter != type_index_lookup.right.end(); ++typeiter) {
530 output << " [" << typeiter->first << "," << typeiter->second << "]";
531 }
532 LOG(3, "DEBUG: type to index map " << output.str());
533 }
534
535 /// 1. construct boost::graph from arguments_t (iter)
536 const size_t N = type_index_lookup.left.size();
537 UndirectedGraph fragmentgraph(N);
538 for(FunctionModel::arguments_t::const_iterator iter = args.begin();
539 iter != args.end(); ++iter) {
540 if (iter->bonded)
541 boost::add_edge(iter->indices.first, iter->indices.second, fragmentgraph);
542 }
543 index_map_t index_map = boost::get(boost::vertex_index, fragmentgraph);
544 LOG(2, "DEBUG: We have " << boost::num_vertices(fragmentgraph) << " nodes in the fragment graph.");
545
546 /// 2. grab first node of the binding model (gives a type)
547 const BindingModel::vector_nodes_t::const_iterator firstiter = _bindingmodel.getNodes().begin();
548 const FragmentNode &firstnode = *firstiter;
549 const Extractors::ParticleType_t &firsttype = firstnode.getAtomicNumber();
550
551 /// 3. grab all candidate nodes contained in arguments_t
552 std::pair<
553 type_index_lookup_t::right_const_iterator,
554 type_index_lookup_t::right_const_iterator> range = type_index_lookup.right.equal_range(firsttype);
555
556 /// 4. go over all candidate nodes (gives index)
557 const size_t nodes_in_bindingmodel = _bindingmodel.getNodes().size();
558 LOG(2, "DEBUG: There are " << nodes_in_bindingmodel << " nodes in the binding model.");
559 set_of_nodes_t set_of_nodes;
560 for (type_index_lookup_t::right_const_iterator rangeiter = range.first;
561 rangeiter != range.second; ++rangeiter) {
562 const size_t &rootindex = rangeiter->second;
563 LOG(2, "DEBUG: Current root index is " << rootindex);
564
565 /// 5a. from node in graph (with this index) perform BFS till n-1 (#nodes in binding model)
566 std::vector<size_t> distances(boost::num_vertices(fragmentgraph));
567 boost::breadth_first_search(
568 fragmentgraph,
569 boost::vertex(rootindex, fragmentgraph),
570 boost::visitor(record_distance(&distances[0])));
571 LOG(3, "DEBUG: BFS discovered the following distances " << distances);
572
573 /// 5b. and store all node in map with distance to root as key
574 nodes_per_level_t nodes_per_level;
575 for (size_t i=0;i<distances.size();++i) {
576 nodes_per_level.insert( std::make_pair(level_t(distances[i]), node_t(i)) );
577 }
578 LOG(3, "DEBUG: Nodes per level are " << nodes_per_level);
579
580 /// 6. construct all possible induced connected subgraphs with this map (see fragmentation)
581 nodes_t nodes;
582 // rootindex is always contained
583 nodes.insert( rootindex );
584 level_t level = 0;
585
586 generateAllInducedConnectedSubgraphs(
587 nodes_in_bindingmodel-1, level, nodes, set_of_nodes, nodes_per_level, fragmentgraph, distances, index_map);
588 LOG(2, "DEBUG: We have found " << set_of_nodes.size() << " unique induced, connected subgraphs.");
589 }
590
591 /// 8. go through each induced, connected subgraph
592 for (set_of_nodes_t::const_iterator nodesiter = set_of_nodes.begin();
593 nodesiter != set_of_nodes.end(); ++nodesiter) {
594 /// 9. convert its nodes into a HomologyGraph using the subgraph and arguments_t (types for index)
595 const nodes_t &nodes = *nodesiter;
596 const HomologyGraph nodes_graph =
597 createHomologyGraphFromNodes(nodes, type_index_lookup, fragmentgraph, index_map);
598
599 /// 10. compare each converted HomologyGraph with _bindingmodel: Accept or Reject
600 if (nodes_graph == _bindingmodel.getGraph()) {
601 LOG(2, "ACCEPT: " << nodes_graph << " is identical to " << _bindingmodel.getGraph());
602 /// 11. for each accepted keyset, pick _all_ symmetric distances from arguments_t
603 FunctionModel::arguments_t argumentbunch;
604 argumentbunch.reserve(args.size());
605 for(FunctionModel::arguments_t::const_iterator iter = args.begin();
606 iter != args.end(); ++iter) {
607 if ((nodes.count(iter->indices.first) != 0) && (nodes.count(iter->indices.second) != 0)) {
608 argumentbunch.push_back(*iter);
609 }
610 }
611 const size_t num_symmetric_distances = nodes.size()*(nodes.size()-1)/2;
612 ASSERT( argumentbunch.size() == num_symmetric_distances,
613 "Extractors::reorderArgumentsByParticleTypes() - only "+toString(argumentbunch.size())+
614 " found instead of "+toString(num_symmetric_distances));
615 LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
616
617 /**
618 * We still need to bring the arguments in the correct order for the potential.
619 *
620 * The potential gives the desired order via the nodes in the HomologyGraph! Their
621 * sequence matches with the user-defined particle types and because of the properties
622 * of the homology graph (FragmentNode stores type and number of edges) it also
623 * matches with the desired bonding graph.
624 */
625
626 /// So, we need to extract from the argumentbunch the information contained in each
627 /// FragmentNode, namely its type and the number of connected edges
628 node_FragmentNode_map_t node_FragmentNode_map = fillNodeFragmentMap(argumentbunch);
629
630 /// Then, we step through the nodes of the bindingmodel ...
631 /// ... and find a suitable mapping from indices in the arguments to
632 /// the index in the order of the HomologyGraph's nodes
633 const BindingModel::vector_nodes_t bindingmodel_nodes = _bindingmodel.getNodes();
634 argindex_to_nodeindex_t argindex_to_nodeindex;
635 size_t nodeindex = 0;
636 for (BindingModel::vector_nodes_t::const_iterator nodeiter = bindingmodel_nodes.begin();
637 nodeiter != bindingmodel_nodes.end(); ++nodeiter) {
638 const FragmentNode &node = *nodeiter;
639 LOG(3, "DEBUG: Binding model's node #" << node << ".");
640 /// ... and pick for each the first (and unique) from these stored nodes.
641 node_FragmentNode_map_t::iterator mapiter = node_FragmentNode_map.begin();
642 for (;mapiter != node_FragmentNode_map.end(); ++mapiter) {
643 if ((mapiter->second.first == node.getAtomicNumber())
644 && (mapiter->second.second == node.getConnectedEdges())) {
645 LOG(3, "DEBUG: #" << mapiter->first << " with type " << mapiter->second.first
646 << " and " << mapiter->second.second << " connected edges matches.");
647 break;
648 }
649 }
650 ASSERT( mapiter != node_FragmentNode_map.end(),
651 "Extractors::reorderArgumentsByParticleTypes() - could not find a suitable node for #"+
652 toString(mapiter->first)+" with type "+toString(mapiter->second.first)+" and "+
653 toString(mapiter->second.second)+" connected edges");
654 std::pair<argindex_to_nodeindex_t::iterator, bool> inserter =
655 argindex_to_nodeindex.insert( std::make_pair(mapiter->first, nodeindex++) );
656 ASSERT( inserter.second,
657 "Extractors::reorderArgumentsByParticleTypes() - node #"+toString(mapiter->first)+
658 " is already present?");
659 // remove to ensure uniqueness of choice
660 node_FragmentNode_map.erase(mapiter);
661 }
662 LOG(4, "DEBUG: argument's indices to node index map is " << argindex_to_nodeindex);
663 // i.e. this is not the arg's index in argumentbunch, but the index of the position
664 // contained in the argument_t
665
666 /// Finally, we only need to bring the arguments in the typical order:
667 /// 01 02 03 04 ... 0n, 12 13 14 ... 1n, 23 24 25 ... 2n, ..., (n-1)n
668 /// These ordering we store in a map for each argument's indices
669 const size_t num_args = argindex_to_nodeindex.size();
670 argument_placement_map_t argument_placement_map = fillArgumentsPlacementMap(num_args);
671 LOG(4, "DEBUG: Placement map is " << argument_placement_map);
672 ASSERT( argument_placement_map.size() == argumentbunch.size(),
673 "Extractors::reorderArgumentsByParticleTypes() - placement map has size "+
674 toString(argument_placement_map.size())+" and we expected "+toString(argumentbunch.size()));
675
676 // and finally resort the arguments with the constructed placement map
677 FunctionModel::arguments_t sortedargs(argumentbunch.size());
678 for (FunctionModel::arguments_t::const_iterator argiter = argumentbunch.begin();
679 argiter != argumentbunch.end(); ++argiter) {
680 const argument_t &arg = *argiter;
681 const argument_t::indices_t translatedIndices =
682 translateIndices(argindex_to_nodeindex, arg.indices);
683 const argument_placement_map_t::const_iterator indexiter =
684 argument_placement_map.find( translatedIndices );
685 ASSERT( indexiter != argument_placement_map.end(),
686 "Extractors::reorderArgumentsByParticleTypes() - could not find place for edge "+
687 toString(arg.indices));
688 sortedargs[indexiter->second] = arg;
689 LOG(3, "DEBUG: Placed argument " << arg << " at #" << indexiter->second
690 << " as translated indices are " << translatedIndices);
691 }
692 LOG(2, "DEBUG: Sorted arguments are " << sortedargs << ".");
693
694 returnargs.push_back(sortedargs);
695 } else {
696 LOG(2, "REJECT: " << nodes_graph << " is not identical to " << _bindingmodel.getGraph());
697 }
698 }
699
700 return returnargs;
701}
702
703FunctionModel::list_of_arguments_t Extractors::filterArgumentsByParticleTypes(
704 const FunctionModel::arguments_t &args,
705 const HomologyGraph &_graph,
706 const ParticleTypes_t &_types,
707 const BindingModel &_bindingmodel
708 )
709{
710 // list allows for quicker removal than vector
711 typedef std::list< argument_t > ListArguments_t;
712 ListArguments_t availableList(args.begin(), args.end());
713 LOG(2, "DEBUG: Initial list of args is " << args << ".");
714
715 // TODO: fill a lookup map such that we don't have O(M^3) scaling, if M is number
716 // of types (and we always must have M(M-1)/2 args) but O(M^2 log(M)). However, as
717 // M is very small (<=3), this is not necessary fruitful now.
718// typedef ParticleTypes_t firsttype;
719// typedef ParticleTypes_t secondtype;
720// typedef std::map< firsttype, std::map< secondtype, FunctionModel::arguments_t > > ArgsLookup_t;
721// ArgsLookup_t ArgsLookup;
722
723 ASSERT( _types.size() <= 2,
724 "Extractors::filterArgumentsByParticleTypes() - this only filters and does not reorder."
725 +std::string(" Hence, it is not useful for multiple arguments per model."));
726
727 // basically, we have two choose any two pairs out of types but only those
728 // where the first is less than the latter. Hence, we start the second
729 // iterator at the current position of the first one and skip the equal case.
730 FunctionModel::list_of_arguments_t returnargs;
731 for (ParticleTypes_t::const_iterator firstiter = _types.begin();
732 firstiter != _types.end();
733 ++firstiter) {
734 for (ParticleTypes_t::const_iterator seconditer = firstiter;
735 seconditer != _types.end();
736 ++seconditer) {
737 if (seconditer == firstiter)
738 continue;
739 LOG(3, "DEBUG: Looking for (" << *firstiter << "," << *seconditer << ") in all args.");
740
741 // search the right one in _args (we might allow switching places of
742 // firstiter and seconditer, as distance is symmetric).
743 ListArguments_t::iterator iter = availableList.begin();
744 while (iter != availableList.end()) {
745 LOG(4, "DEBUG: Current args is " << *iter << ".");
746 if ((iter->types.first == *firstiter)
747 && (iter->types.second == *seconditer)) {
748 returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
749 iter = availableList.erase(iter);
750 LOG(4, "DEBUG: Accepted argument.");
751 } else if ((iter->types.first == *seconditer)
752 && (iter->types.second == *firstiter)) {
753 returnargs.push_back( FunctionModel::arguments_t(1, *iter) );
754 iter = availableList.erase(iter);
755 LOG(4, "DEBUG: Accepted (flipped) argument.");
756 } else {
757 ++iter;
758 LOG(4, "DEBUG: Rejected argument.");
759 }
760 }
761 }
762 }
763
764 LOG(2, "DEBUG: We have generated " << returnargs.size() << " tuples of distances.");
765
766 return returnargs;
767}
768
769
770FunctionModel::arguments_t Extractors::combineArguments(
771 const FunctionModel::arguments_t &firstargs,
772 const FunctionModel::arguments_t &secondargs)
773{
774 FunctionModel::arguments_t args = concatenateArguments(firstargs, secondargs);
775 std::sort(args.begin(), args.end(),
776 boost::bind(&argument_t::operator<, _1, _2));
777 FunctionModel::arguments_t::iterator iter =
778 std::unique(args.begin(), args.end(),
779 boost::bind(&argument_t::operator==, _1, _2));
780 args.erase(iter, args.end());
781 return args;
782}
783
784FunctionModel::arguments_t Extractors::concatenateArguments(
785 const FunctionModel::arguments_t &firstargs,
786 const FunctionModel::arguments_t &secondargs)
787{
788 FunctionModel::arguments_t args(firstargs);
789 args.insert(args.end(), secondargs.begin(), secondargs.end());
790 return args;
791}
792
793FunctionModel::list_of_arguments_t Extractors::concatenateListOfArguments(
794 const FunctionModel::list_of_arguments_t &firstlistargs,
795 const FunctionModel::list_of_arguments_t &secondlistargs)
796{
797 FunctionModel::list_of_arguments_t listargs(firstlistargs);
798 listargs.insert(listargs.end(), secondlistargs.begin(), secondlistargs.end());
799 return listargs;
800}
Note: See TracBrowser for help on using the repository browser.