source: src/FunctionApproximation/Extractors.cpp@ 0b7036

Fix_FitPotential_needs_atomicnumbers
Last change on this file since 0b7036 was 0b7036, checked in by Frederik Heber <heber@…>, 9 years ago

Added getSubgraph() to all EmpiricalPotentials.

  • this encodes the specific bonding structure which this potential represents and which allows finding all(!) subgraphs within the fragment graph in the correct manner.
  • Property mode set to 100644
File size: 16.8 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 <sstream>
40#include <utility>
41#include <vector>
42#include <boost/assign.hpp>
43#include <boost/bind.hpp>
44#include <boost/foreach.hpp>
45#include <boost/graph/adjacency_list.hpp>
46
47#include "CodePatterns/Assert.hpp"
48#include "CodePatterns/IteratorAdaptors.hpp"
49#include "CodePatterns/Log.hpp"
50#include "CodePatterns/toString.hpp"
51
52#include "LinearAlgebra/Vector.hpp"
53
54#include "FunctionApproximation/Extractors.hpp"
55#include "FunctionApproximation/FunctionArgument.hpp"
56
57using namespace boost::assign;
58
59FunctionModel::arguments_t
60Extractors::gatherAllSymmetricDistanceArguments(
61 const Fragment::positions_t& positions,
62 const Fragment::atomicnumbers_t& atomicnumbers,
63 const size_t globalid)
64{
65 FunctionModel::arguments_t result;
66
67 // go through current configuration and gather all other distances
68 Fragment::positions_t::const_iterator firstpositer = positions.begin();
69 for (;firstpositer != positions.end(); ++firstpositer) {
70 Fragment::positions_t::const_iterator secondpositer = firstpositer;
71 for (; secondpositer != positions.end(); ++secondpositer) {
72 if (firstpositer == secondpositer)
73 continue;
74 argument_t arg;
75 const Vector firsttemp((*firstpositer)[0],(*firstpositer)[1],(*firstpositer)[2]);
76 const Vector secondtemp((*secondpositer)[0],(*secondpositer)[1],(*secondpositer)[2]);
77 arg.distance = firsttemp.distance(secondtemp);
78 arg.types = std::make_pair(
79 (int)atomicnumbers[ std::distance(positions.begin(), firstpositer) ],
80 (int)atomicnumbers[ std::distance(positions.begin(), secondpositer) ]
81 );
82 arg.indices = std::make_pair(
83 std::distance(
84 positions.begin(), firstpositer),
85 std::distance(
86 positions.begin(), secondpositer)
87 );
88 arg.globalid = globalid;
89 LOG(3, "DEBUG: Created argument " << arg << ".");
90 result.push_back(arg);
91 }
92 }
93
94 return result;
95}
96
97Extractors::elementcounts_t
98Extractors::_detail::getElementCounts(
99 const Fragment::atomicnumbers_t elements
100 )
101{
102 elementcounts_t elementcounts;
103 for (Fragment::atomicnumbers_t::const_iterator elementiter = elements.begin();
104 elementiter != elements.end(); ++elementiter) {
105 // insert new element
106 std::pair< elementcounts_t::iterator, bool> inserter =
107 elementcounts.insert( std::make_pair( *elementiter, 1) );
108 // if already present, just increase its count
109 if (!inserter.second)
110 ++(inserter.first->second);
111 }
112 return elementcounts;
113}
114
115struct ParticleTypesComparator {
116 bool operator()(const argument_t::types_t &a, const argument_t::types_t &b)
117 {
118 if (a.first < a.second) {
119 if (b.first < b.second) {
120 if (a.first < b.first)
121 return true;
122 else if (a.first > b.first)
123 return false;
124 else
125 return (a.second < b.second);
126 } else {
127 if (a.first < b.second)
128 return true;
129 else if (a.first > b.second)
130 return false;
131 else
132 return (a.second < b.first);
133 }
134 } else {
135 if (b.first < b.second) {
136 if (a.second < b.first)
137 return true;
138 else if (a.second > b.first)
139 return false;
140 else
141 return (a.first < b.second);
142 } else {
143 if (a.second < b.second)
144 return true;
145 else if (a.second > b.second)
146 return false;
147 else
148 return (a.first < b.first);
149 }
150 }
151 }
152};
153
154std::ostream& operator<<(std::ostream &out, const argument_t::types_t &a)
155{
156 out << "[" << a.first << "," << a.second << "]";
157 return out;
158}
159
160FunctionModel::list_of_arguments_t Extractors::reorderArgumentsByParticleTypes(
161 const FunctionModel::list_of_arguments_t &listargs,
162 const ParticleTypes_t &_types
163 )
164{
165 FunctionModel::list_of_arguments_t returnargs;
166 for (FunctionModel::list_of_arguments_t::const_iterator iter = listargs.begin();
167 iter != listargs.end(); ++iter) {
168 const FunctionModel::arguments_t &args = *iter;
169 /// We place all arguments into multimap according to particle type pair.
170 // here, we need a special comparator such that types in key pair are always
171 // properly ordered.
172 typedef std::multimap<
173 argument_t::types_t,
174 argument_t,
175 ParticleTypesComparator> TypePair_Argument_Map_t;
176 TypePair_Argument_Map_t argument_map;
177 for(FunctionModel::arguments_t::const_iterator iter = args.begin();
178 iter != args.end(); ++iter) {
179 argument_map.insert( std::make_pair(iter->types, *iter) );
180 }
181 LOG(4, "DEBUG: particle_type map is " << argument_map << ".");
182
183 /// Then, we create the desired unique keys
184 typedef std::vector<argument_t::types_t> UniqueTypes_t;
185 UniqueTypes_t UniqueTypes;
186 for (ParticleTypes_t::const_iterator firstiter = _types.begin();
187 firstiter != _types.end();
188 ++firstiter) {
189 for (ParticleTypes_t::const_iterator seconditer = firstiter;
190 seconditer != _types.end();
191 ++seconditer) {
192 if (seconditer == firstiter)
193 continue;
194 UniqueTypes.push_back( std::make_pair(*firstiter, *seconditer) );
195 }
196 }
197 LOG(4, "DEBUG: Created unique types as keys " << UniqueTypes << ".");
198
199 /// turn fragment into graph from the known edges
200 boost::adjacency_list<> fragmentgraph(_types.size());
201 for(FunctionModel::arguments_t::const_iterator iter = args.begin();
202 iter != args.end(); ++iter)
203 boost::add_edge(iter->indices.first, iter->indices.second, fragmentgraph);
204 boost::graph_traits < boost::adjacency_list <> >::vertex_iterator i, end;
205 boost::graph_traits < boost::adjacency_list <> >::adjacency_iterator ai, a_end;
206 boost::property_map < boost::adjacency_list <>, boost::vertex_index_t >::type
207 index_map = boost::get(boost::vertex_index, fragmentgraph);
208
209 /// Finally, we use the unique key list to pick corresponding arguments from the map
210 FunctionModel::arguments_t sortedargs;
211 sortedargs.reserve(args.size());
212 while (!argument_map.empty()) {
213 // note that particle_types_t may be flipped, i.e. 1,8 is equal to 8,1, but we
214 // must maintain the correct order in indices in accordance with the order
215 // in _types, i.e. 1,8,1 must match with e.g. ids 1,0,2 where 1 has type 1,
216 // 0 has type 8, and 2 has type 2.
217 // In other words: We do not want to flip/modify arguments such that they match
218 // with the specific type pair we seek but then this comes at the price that we
219 // have flip indices when the types in a pair are flipped.
220
221 typedef std::vector<size_t> indices_t;
222 //!> here, we gather the indices as we discover them
223 indices_t indices;
224 indices.resize(_types.size(), (size_t)-1);
225
226 // these are two iterators that create index pairs in the same way as we have
227 // created type pairs. If a -1 is still present in indices, then the index is
228 // still arbitrary but is then set by the next found index
229 indices_t::iterator firstindex = indices.begin();
230 indices_t::iterator secondindex = firstindex+1;
231
232 //!> here, we gather the current bunch of arguments as we find them
233 FunctionModel::arguments_t argumentbunch;
234 argumentbunch.reserve(UniqueTypes.size());
235
236 for (UniqueTypes_t::const_iterator typeiter = UniqueTypes.begin();
237 typeiter != UniqueTypes.end(); ++typeiter) {
238 // have all arguments to same type pair as list within the found range
239 std::pair<
240 TypePair_Argument_Map_t::iterator,
241 TypePair_Argument_Map_t::iterator> range_t =
242 argument_map.equal_range(*typeiter);
243 LOG(4, "DEBUG: Set of arguments to current key [" << typeiter->first << ","
244 << typeiter->second << "] is " << std::list<argument_t>(
245 MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.first),
246 MapValueIterator<TypePair_Argument_Map_t::iterator>(range_t.second)
247 ) << ".");
248 // the first key is always easy and is pivot which the rest has to be associated to
249 if (typeiter == UniqueTypes.begin()) {
250 const argument_t & arg = range_t.first->second;
251 if ((typeiter->first == arg.types.first) && (typeiter->second == arg.types.second)) {
252 // store in correct order
253 *firstindex = arg.indices.first;
254 *secondindex = arg.indices.second;
255 } else {
256 // store in flipped order
257 *firstindex = arg.indices.second;
258 *secondindex = arg.indices.first;
259 }
260 argumentbunch.push_back(arg);
261 argument_map.erase(range_t.first);
262 LOG(4, "DEBUG: Gathered first argument " << arg << ".");
263 } else {
264 // go through the range and pick the first argument matching the index constraints
265 for (TypePair_Argument_Map_t::iterator argiter = range_t.first;
266 argiter != range_t.second; ++argiter) {
267 // seconditer may be -1 still
268 const argument_t &arg = argiter->second;
269 if (arg.indices.first == *firstindex) {
270 if ((arg.indices.second == *secondindex) || (*secondindex == (size_t)-1)) {
271 if (*secondindex == (size_t)-1)
272 *secondindex = arg.indices.second;
273 argumentbunch.push_back(arg);
274 argument_map.erase(argiter);
275 LOG(4, "DEBUG: Gathered another argument " << arg << ".");
276 break;
277 }
278 } else if ((arg.indices.first == *secondindex) || (*secondindex == (size_t)-1)) {
279 if (arg.indices.second == *firstindex) {
280 if (*secondindex == (size_t)-1)
281 *secondindex = arg.indices.first;
282 argumentbunch.push_back(arg);
283 argument_map.erase(argiter);
284 LOG(4, "DEBUG: Gathered another (flipped) argument " << arg << ".");
285 break;
286 }
287 }
288 }
289 }
290 // move along in indices and check bounds
291 ++secondindex;
292 if (secondindex == indices.end()) {
293 ++firstindex;
294 if (firstindex != indices.end()-1)
295 secondindex = firstindex+1;
296 }
297 }
298 ASSERT( (firstindex == indices.end()-1) && (secondindex == indices.end()),
299 "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
300 ASSERT( argumentbunch.size() == UniqueTypes.size(),
301 "Extractors::reorderArgumentsByParticleTypes() - we have not gathered enough arguments.");
302 // place bunch of arguments in return args
303 LOG(3, "DEBUG: Given types " << _types << " and found indices " << indices << ".");
304 LOG(3, "DEBUG: Final bunch of arguments is " << argumentbunch << ".");
305 sortedargs.insert(sortedargs.end(), argumentbunch.begin(), argumentbunch.end());
306 }
307 returnargs.push_back(sortedargs);
308 }
309
310 return returnargs;
311}
312
313FunctionModel::list_of_arguments_t Extractors::filterArgumentsByParticleTypes(
314 const FunctionModel::arguments_t &args,
315 const ParticleTypes_t &_types
316 )
317{
318 typedef std::list< argument_t > ListArguments_t;
319 ListArguments_t availableList(args.begin(), args.end());
320 LOG(2, "DEBUG: Initial list of args is " << args << ".");
321
322
323 // TODO: fill a lookup map such that we don't have O(M^3) scaling, if M is number
324 // of types (and we always must have M(M-1)/2 args) but O(M^2 log(M)). However, as
325 // M is very small (<=3), this is not necessary fruitful now.
326// typedef ParticleTypes_t firsttype;
327// typedef ParticleTypes_t secondtype;
328// typedef std::map< firsttype, std::map< secondtype, boost::ref(args) > > ArgsLookup_t;
329// ArgsLookup_t ArgsLookup;
330
331 // basically, we have two choose any two pairs out of types but only those
332 // where the first is less than the latter. Hence, we start the second
333 // iterator at the current position of the first one and skip the equal case.
334 FunctionModel::arguments_t allargs;
335 allargs.reserve(args.size());
336 for (ParticleTypes_t::const_iterator firstiter = _types.begin();
337 firstiter != _types.end();
338 ++firstiter) {
339 for (ParticleTypes_t::const_iterator seconditer = firstiter;
340 seconditer != _types.end();
341 ++seconditer) {
342 if (seconditer == firstiter)
343 continue;
344 LOG(3, "DEBUG: Looking for (" << *firstiter << "," << *seconditer << ") in all args.");
345
346 // search the right one in _args (we might allow switching places of
347 // firstiter and seconditer, as distance is symmetric).
348 ListArguments_t::iterator iter = availableList.begin();
349 while (iter != availableList.end()) {
350 LOG(4, "DEBUG: Current args is " << *iter << ".");
351 if ((iter->types.first == *firstiter)
352 && (iter->types.second == *seconditer)) {
353 allargs.push_back( *iter );
354 iter = availableList.erase(iter);
355 LOG(4, "DEBUG: Accepted argument.");
356 } else if ((iter->types.first == *seconditer)
357 && (iter->types.second == *firstiter)) {
358 allargs.push_back( *iter );
359 iter = availableList.erase(iter);
360 LOG(4, "DEBUG: Accepted (flipped) argument.");
361 } else {
362 ++iter;
363 LOG(4, "DEBUG: Rejected argument.");
364 }
365 }
366 }
367 }
368 LOG(2, "DEBUG: Final list of args is " << allargs << ".");
369
370 // first, we bring together tuples of distances that belong together
371 FunctionModel::list_of_arguments_t singlelist_allargs;
372 singlelist_allargs.push_back(allargs);
373 FunctionModel::list_of_arguments_t sortedargs =
374 reorderArgumentsByParticleTypes(singlelist_allargs, _types);
375 ASSERT( sortedargs.size() == (size_t)1,
376 "Extractors::filterArgumentsByParticleTypes() - reordering did not generate a single list.");
377 // then we split up the tuples of arguments and place each into single list
378 FunctionModel::list_of_arguments_t returnargs;
379 FunctionModel::arguments_t::const_iterator argiter = sortedargs.begin()->begin();
380 const size_t num_types = _types.size();
381 const size_t args_per_tuple = num_types * (num_types-1) / 2;
382 while (argiter != sortedargs.begin()->end()) {
383 FunctionModel::arguments_t currenttuple(args_per_tuple);
384 const FunctionModel::arguments_t::const_iterator startiter = argiter;
385 std::advance(argiter, args_per_tuple);
386#ifndef NDEBUG
387 FunctionModel::arguments_t::const_iterator endoutiter =
388#endif
389 std::copy(startiter, argiter, currenttuple.begin());
390 ASSERT( endoutiter == currenttuple.end(),
391 "Extractors::filterArgumentsByParticleTypes() - currenttuple not initialized to right size.");
392 returnargs.push_back(currenttuple);
393 }
394
395 LOG(2, "DEBUG: We have generated " << returnargs.size() << " tuples of distances.");
396
397 return returnargs;
398}
399
400
401FunctionModel::arguments_t Extractors::combineArguments(
402 const FunctionModel::arguments_t &firstargs,
403 const FunctionModel::arguments_t &secondargs)
404{
405 FunctionModel::arguments_t args = concatenateArguments(firstargs, secondargs);
406 std::sort(args.begin(), args.end(),
407 boost::bind(&argument_t::operator<, _1, _2));
408 FunctionModel::arguments_t::iterator iter =
409 std::unique(args.begin(), args.end(),
410 boost::bind(&argument_t::operator==, _1, _2));
411 args.erase(iter, args.end());
412 return args;
413}
414
415FunctionModel::arguments_t Extractors::concatenateArguments(
416 const FunctionModel::arguments_t &firstargs,
417 const FunctionModel::arguments_t &secondargs)
418{
419 FunctionModel::arguments_t args(firstargs);
420 args.insert(args.end(), secondargs.begin(), secondargs.end());
421 return args;
422}
423
424FunctionModel::list_of_arguments_t Extractors::concatenateListOfArguments(
425 const FunctionModel::list_of_arguments_t &firstlistargs,
426 const FunctionModel::list_of_arguments_t &secondlistargs)
427{
428 FunctionModel::list_of_arguments_t listargs(firstlistargs);
429 listargs.insert(listargs.end(), secondlistargs.begin(), secondlistargs.end());
430 return listargs;
431}
Note: See TracBrowser for help on using the repository browser.