source: src/FunctionApproximation/Extractors.cpp@ ffad62

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

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