1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team.
5 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
6 * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS 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 GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37 /*! \internal \file
38 * \brief
39 * Implements functions in indexutil.h.
40 *
41 * \author Teemu Murtola <teemu.murtola@gmail.com>
42 * \ingroup module_selection
43 */
44 #include "gmxpre.h"
45
46 #include "indexutil.h"
47
48 #include <cstdlib>
49 #include <cstring>
50
51 #include <algorithm>
52 #include <numeric>
53 #include <string>
54 #include <vector>
55
56 #include "gromacs/topology/block.h"
57 #include "gromacs/topology/index.h"
58 #include "gromacs/topology/mtop_lookup.h"
59 #include "gromacs/topology/mtop_util.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/utility/exceptions.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/smalloc.h"
64 #include "gromacs/utility/stringutil.h"
65 #include "gromacs/utility/textwriter.h"
66
67 namespace gmx
68 {
69
IndexGroupsAndNames(const t_blocka & indexGroup,ArrayRef<char const * const> groupNames)70 IndexGroupsAndNames::IndexGroupsAndNames(const t_blocka& indexGroup, ArrayRef<char const* const> groupNames) :
71 indexGroup_{ indexGroup }
72 {
73 std::copy(groupNames.begin(), groupNames.end(), std::back_inserter(groupNames_));
74 GMX_ASSERT(indexGroup_.nr == ssize(groupNames),
75 "Number of groups must match number of group names.");
76 }
77
containsGroupName(const std::string & groupName) const78 bool IndexGroupsAndNames::containsGroupName(const std::string& groupName) const
79 {
80 return std::any_of(
81 std::begin(groupNames_), std::end(groupNames_),
82 [&groupName](const std::string& name) { return equalCaseInsensitive(groupName, name); });
83 }
84
indices(const std::string & groupName) const85 std::vector<index> IndexGroupsAndNames::indices(const std::string& groupName) const
86 {
87 if (!containsGroupName(groupName))
88 {
89 GMX_THROW(
90 InconsistentInputError(
91 std::string("Group ") + groupName
92 + " referenced in the .mdp file was not found in the index file.\n"
93 "Group names must match either [moleculetype] names or custom index "
94 "group\n"
95 "names, in which case you must supply an index file to the '-n' option\n"
96 "of grompp."));
97 }
98 const auto groupNamePosition = std::find_if(
99 std::begin(groupNames_), std::end(groupNames_),
100 [&groupName](const std::string& name) { return equalCaseInsensitive(groupName, name); });
101 const auto groupIndex = std::distance(std::begin(groupNames_), groupNamePosition);
102 const auto groupSize = indexGroup_.index[groupIndex + 1] - indexGroup_.index[groupIndex];
103 std::vector<index> groupIndices(groupSize);
104 const auto startingIndex = indexGroup_.index[groupIndex];
105 std::iota(std::begin(groupIndices), std::end(groupIndices), startingIndex);
106 std::transform(std::begin(groupIndices), std::end(groupIndices), std::begin(groupIndices),
107 [blockLookup = indexGroup_.a](auto i) { return blockLookup[i]; });
108 return groupIndices;
109 }
110
111 } // namespace gmx
112
113 /********************************************************************
114 * gmx_ana_indexgrps_t functions
115 ********************************************************************/
116
117 /*! \internal \brief
118 * Stores a set of index groups.
119 */
120 struct gmx_ana_indexgrps_t
121 {
122 //! Initializes an empty set of groups.
gmx_ana_indexgrps_tgmx_ana_indexgrps_t123 explicit gmx_ana_indexgrps_t(int nr) : g(nr) { names.reserve(nr); }
~gmx_ana_indexgrps_tgmx_ana_indexgrps_t124 ~gmx_ana_indexgrps_t()
125 {
126 for (auto& indexGrp : g)
127 {
128 gmx_ana_index_deinit(&indexGrp);
129 }
130 }
131
132
133 /** Array of index groups. */
134 std::vector<gmx_ana_index_t> g;
135 /** Group names. */
136 std::vector<std::string> names;
137 };
138
139 /*!
140 * \param[out] g Index group structure.
141 * \param[in] top Topology structure.
142 * \param[in] fnm File name for the index file.
143 * Memory is automatically allocated.
144 *
145 * One or both of \p top or \p fnm can be NULL.
146 * If \p top is NULL, an index file is required and the groups are read
147 * from the file (uses Gromacs routine init_index()).
148 * If \p fnm is NULL, default groups are constructed based on the
149 * topology (uses Gromacs routine analyse()).
150 * If both are null, the index group structure is initialized empty.
151 */
gmx_ana_indexgrps_init(gmx_ana_indexgrps_t ** g,gmx_mtop_t * top,const char * fnm)152 void gmx_ana_indexgrps_init(gmx_ana_indexgrps_t** g, gmx_mtop_t* top, const char* fnm)
153 {
154 t_blocka* block = nullptr;
155 char** names = nullptr;
156
157 if (fnm)
158 {
159 block = init_index(fnm, &names);
160 }
161 else if (top)
162 {
163 block = new_blocka();
164 // TODO: Propagate mtop further.
165 t_atoms atoms = gmx_mtop_global_atoms(top);
166 analyse(&atoms, block, &names, FALSE, FALSE);
167 done_atom(&atoms);
168 }
169 else
170 {
171 *g = new gmx_ana_indexgrps_t(0);
172 return;
173 }
174
175 try
176 {
177 *g = new gmx_ana_indexgrps_t(block->nr);
178 for (int i = 0; i < block->nr; ++i)
179 {
180 gmx_ana_index_t* grp = &(*g)->g[i];
181
182 grp->isize = block->index[i + 1] - block->index[i];
183 snew(grp->index, grp->isize);
184 for (int j = 0; j < grp->isize; ++j)
185 {
186 grp->index[j] = block->a[block->index[i] + j];
187 }
188 grp->nalloc_index = grp->isize;
189 (*g)->names.emplace_back(names[i]);
190 }
191 }
192 catch (...)
193 {
194 for (int i = 0; i < block->nr; ++i)
195 {
196 sfree(names[i]);
197 }
198 sfree(names);
199 done_blocka(block);
200 sfree(block);
201 throw;
202 }
203 for (int i = 0; i < block->nr; ++i)
204 {
205 sfree(names[i]);
206 }
207 sfree(names);
208 done_blocka(block);
209 sfree(block);
210 }
211
212 /*!
213 * \param[in] g Index groups structure.
214 *
215 * The pointer \p g is invalid after the call.
216 */
gmx_ana_indexgrps_free(gmx_ana_indexgrps_t * g)217 void gmx_ana_indexgrps_free(gmx_ana_indexgrps_t* g)
218 {
219 delete g;
220 }
221
222
223 /*!
224 * \param[out] dest Output structure.
225 * \param[out] destName Receives the name of the group if found.
226 * \param[in] src Input index groups.
227 * \param[in] n Number of the group to extract.
228 * \returns true if \p n is a valid group in \p src, false otherwise.
229 */
gmx_ana_indexgrps_extract(gmx_ana_index_t * dest,std::string * destName,gmx_ana_indexgrps_t * src,int n)230 bool gmx_ana_indexgrps_extract(gmx_ana_index_t* dest, std::string* destName, gmx_ana_indexgrps_t* src, int n)
231 {
232 destName->clear();
233 if (n < 0 || n >= gmx::index(src->g.size()))
234 {
235 dest->isize = 0;
236 return false;
237 }
238
239 if (destName != nullptr)
240 {
241 *destName = src->names[n];
242 }
243 gmx_ana_index_copy(dest, &src->g[n], true);
244 return true;
245 }
246
247 /*!
248 * \param[out] dest Output structure.
249 * \param[out] destName Receives the name of the group if found.
250 * \param[in] src Input index groups.
251 * \param[in] name Name (or part of the name) of the group to extract.
252 * \returns true if \p name is a valid group in \p src, false otherwise.
253 *
254 * Uses the Gromacs routine find_group() to find the actual group;
255 * the comparison is case-insensitive.
256 */
gmx_ana_indexgrps_find(gmx_ana_index_t * dest,std::string * destName,gmx_ana_indexgrps_t * src,const char * name)257 bool gmx_ana_indexgrps_find(gmx_ana_index_t* dest, std::string* destName, gmx_ana_indexgrps_t* src, const char* name)
258 {
259 const char** names;
260
261 destName->clear();
262 snew(names, src->g.size());
263 for (size_t i = 0; i < src->g.size(); ++i)
264 {
265 names[i] = src->names[i].c_str();
266 }
267 int n = find_group(const_cast<char*>(name), src->g.size(), const_cast<char**>(names));
268 sfree(names);
269 if (n < 0)
270 {
271 dest->isize = 0;
272 return false;
273 }
274
275 return gmx_ana_indexgrps_extract(dest, destName, src, n);
276 }
277
278 /*!
279 * \param[in] writer Writer to use for output.
280 * \param[in] g Index groups to print.
281 * \param[in] maxn Maximum number of indices to print
282 * (-1 = print all, 0 = print only names).
283 */
gmx_ana_indexgrps_print(gmx::TextWriter * writer,gmx_ana_indexgrps_t * g,int maxn)284 void gmx_ana_indexgrps_print(gmx::TextWriter* writer, gmx_ana_indexgrps_t* g, int maxn)
285 {
286 for (gmx::index i = 0; i < gmx::ssize(g->g); ++i)
287 {
288 writer->writeString(gmx::formatString(" Group %2zd \"%s\" ", i, g->names[i].c_str()));
289 gmx_ana_index_dump(writer, &g->g[i], maxn);
290 }
291 }
292
293 /********************************************************************
294 * gmx_ana_index_t functions
295 ********************************************************************/
296
297 /*!
298 * \param[in,out] g Index group structure.
299 * \param[in] isize Maximum number of atoms to reserve space for.
300 */
gmx_ana_index_reserve(gmx_ana_index_t * g,int isize)301 void gmx_ana_index_reserve(gmx_ana_index_t* g, int isize)
302 {
303 if (g->nalloc_index < isize)
304 {
305 srenew(g->index, isize);
306 g->nalloc_index = isize;
307 }
308 }
309
310 /*!
311 * \param[in,out] g Index group structure.
312 *
313 * Resizes the memory allocated for holding the indices such that the
314 * current contents fit.
315 */
gmx_ana_index_squeeze(gmx_ana_index_t * g)316 void gmx_ana_index_squeeze(gmx_ana_index_t* g)
317 {
318 srenew(g->index, g->isize);
319 g->nalloc_index = g->isize;
320 }
321
322 /*!
323 * \param[out] g Output structure.
324 *
325 * Any contents of \p g are discarded without freeing.
326 */
gmx_ana_index_clear(gmx_ana_index_t * g)327 void gmx_ana_index_clear(gmx_ana_index_t* g)
328 {
329 g->isize = 0;
330 g->index = nullptr;
331 g->nalloc_index = 0;
332 }
333
334 /*!
335 * \param[out] g Output structure.
336 * \param[in] isize Number of atoms in the new group.
337 * \param[in] index Array of \p isize atoms (can be NULL if \p isize is 0).
338 * \param[in] nalloc Number of elements allocated for \p index
339 * (if 0, \p index is not freed in gmx_ana_index_deinit())
340 *
341 * No copy if \p index is made.
342 */
gmx_ana_index_set(gmx_ana_index_t * g,int isize,int * index,int nalloc)343 void gmx_ana_index_set(gmx_ana_index_t* g, int isize, int* index, int nalloc)
344 {
345 g->isize = isize;
346 g->index = index;
347 g->nalloc_index = nalloc;
348 }
349
350 /*!
351 * \param[out] g Output structure.
352 * \param[in] natoms Number of atoms.
353 */
gmx_ana_index_init_simple(gmx_ana_index_t * g,int natoms)354 void gmx_ana_index_init_simple(gmx_ana_index_t* g, int natoms)
355 {
356 int i;
357
358 g->isize = natoms;
359 snew(g->index, natoms);
360 for (i = 0; i < natoms; ++i)
361 {
362 g->index[i] = i;
363 }
364 g->nalloc_index = natoms;
365 }
366
367 /*!
368 * \param[in] g Index group structure.
369 *
370 * The pointer \p g is not freed.
371 */
gmx_ana_index_deinit(gmx_ana_index_t * g)372 void gmx_ana_index_deinit(gmx_ana_index_t* g)
373 {
374 if (g->nalloc_index > 0)
375 {
376 sfree(g->index);
377 }
378 gmx_ana_index_clear(g);
379 }
380
381 /*!
382 * \param[out] dest Destination index group.
383 * \param[in] src Source index group.
384 * \param[in] bAlloc If true, memory is allocated at \p dest; otherwise,
385 * it is assumed that enough memory has been allocated for index.
386 */
gmx_ana_index_copy(gmx_ana_index_t * dest,gmx_ana_index_t * src,bool bAlloc)387 void gmx_ana_index_copy(gmx_ana_index_t* dest, gmx_ana_index_t* src, bool bAlloc)
388 {
389 dest->isize = src->isize;
390 if (bAlloc)
391 {
392 snew(dest->index, dest->isize);
393 dest->nalloc_index = dest->isize;
394 }
395 if (dest->isize > 0)
396 {
397 std::memcpy(dest->index, src->index, dest->isize * sizeof(*dest->index));
398 }
399 }
400
401 /*!
402 * \param[in] writer Writer to use for output.
403 * \param[in] g Index group to print.
404 * \param[in] maxn Maximum number of indices to print (-1 = print all).
405 */
gmx_ana_index_dump(gmx::TextWriter * writer,gmx_ana_index_t * g,int maxn)406 void gmx_ana_index_dump(gmx::TextWriter* writer, gmx_ana_index_t* g, int maxn)
407 {
408 writer->writeString(gmx::formatString("(%d atoms)", g->isize));
409 if (maxn != 0)
410 {
411 writer->writeString(":");
412 int n = g->isize;
413 if (maxn >= 0 && n > maxn)
414 {
415 n = maxn;
416 }
417 for (int j = 0; j < n; ++j)
418 {
419 writer->writeString(gmx::formatString(" %d", g->index[j] + 1));
420 }
421 if (n < g->isize)
422 {
423 writer->writeString(" ...");
424 }
425 }
426 writer->ensureLineBreak();
427 }
428
gmx_ana_index_get_max_index(gmx_ana_index_t * g)429 int gmx_ana_index_get_max_index(gmx_ana_index_t* g)
430 {
431 if (g->isize == 0)
432 {
433 return 0;
434 }
435 else
436 {
437 return *std::max_element(g->index, g->index + g->isize);
438 }
439 }
440
441 /*!
442 * \param[in] g Index group to check.
443 * \returns true if the index group is sorted and has no duplicates,
444 * false otherwise.
445 */
gmx_ana_index_check_sorted(gmx_ana_index_t * g)446 bool gmx_ana_index_check_sorted(gmx_ana_index_t* g)
447 {
448 int i;
449
450 for (i = 0; i < g->isize - 1; ++i)
451 {
452 if (g->index[i + 1] <= g->index[i])
453 {
454 return false;
455 }
456 }
457 return true;
458 }
459
gmx_ana_index_check_range(gmx_ana_index_t * g,int natoms)460 bool gmx_ana_index_check_range(gmx_ana_index_t* g, int natoms)
461 {
462 for (int i = 0; i < g->isize; ++i)
463 {
464 if (g->index[i] < 0 || g->index[i] >= natoms)
465 {
466 return false;
467 }
468 }
469 return true;
470 }
471
472 /********************************************************************
473 * Set operations
474 ********************************************************************/
475
476 /*!
477 * \param[in,out] g Index group to be sorted.
478 */
gmx_ana_index_sort(gmx_ana_index_t * g)479 void gmx_ana_index_sort(gmx_ana_index_t* g)
480 {
481 std::sort(g->index, g->index + g->isize);
482 }
483
gmx_ana_index_remove_duplicates(gmx_ana_index_t * g)484 void gmx_ana_index_remove_duplicates(gmx_ana_index_t* g)
485 {
486 int j = 0;
487 for (int i = 0; i < g->isize; ++i)
488 {
489 if (i == 0 || g->index[i - 1] != g->index[i])
490 {
491 g->index[j] = g->index[i];
492 ++j;
493 }
494 }
495 g->isize = j;
496 }
497
498 /*!
499 * \param[in] a Index group to check.
500 * \param[in] b Index group to check.
501 * \returns true if \p a and \p b are equal, false otherwise.
502 */
gmx_ana_index_equals(gmx_ana_index_t * a,gmx_ana_index_t * b)503 bool gmx_ana_index_equals(gmx_ana_index_t* a, gmx_ana_index_t* b)
504 {
505 int i;
506
507 if (a->isize != b->isize)
508 {
509 return false;
510 }
511 for (i = 0; i < a->isize; ++i)
512 {
513 if (a->index[i] != b->index[i])
514 {
515 return false;
516 }
517 }
518 return true;
519 }
520
521 /*!
522 * \param[in] a Index group to check against.
523 * \param[in] b Index group to check.
524 * \returns true if \p b is contained in \p a,
525 * false otherwise.
526 *
527 * If the elements are not in the same order in both groups, the function
528 * fails. However, the groups do not need to be sorted.
529 */
gmx_ana_index_contains(gmx_ana_index_t * a,gmx_ana_index_t * b)530 bool gmx_ana_index_contains(gmx_ana_index_t* a, gmx_ana_index_t* b)
531 {
532 int i, j;
533
534 for (i = j = 0; j < b->isize; ++i, ++j)
535 {
536 while (i < a->isize && a->index[i] != b->index[j])
537 {
538 ++i;
539 }
540 if (i == a->isize)
541 {
542 return false;
543 }
544 }
545 return true;
546 }
547
548 /*!
549 * \param[out] dest Output index group (the intersection of \p a and \p b).
550 * \param[in] a First index group.
551 * \param[in] b Second index group.
552 *
553 * \p dest can be the same as \p a or \p b.
554 */
gmx_ana_index_intersection(gmx_ana_index_t * dest,gmx_ana_index_t * a,gmx_ana_index_t * b)555 void gmx_ana_index_intersection(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b)
556 {
557 int i, j, k;
558
559 for (i = j = k = 0; i < a->isize && j < b->isize; ++i)
560 {
561 while (j < b->isize && b->index[j] < a->index[i])
562 {
563 ++j;
564 }
565 if (j < b->isize && b->index[j] == a->index[i])
566 {
567 dest->index[k++] = b->index[j++];
568 }
569 }
570 dest->isize = k;
571 }
572
573 /*!
574 * \param[out] dest Output index group (the difference \p a - \p b).
575 * \param[in] a First index group.
576 * \param[in] b Second index group.
577 *
578 * \p dest can equal \p a, but not \p b.
579 */
gmx_ana_index_difference(gmx_ana_index_t * dest,gmx_ana_index_t * a,gmx_ana_index_t * b)580 void gmx_ana_index_difference(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b)
581 {
582 int i, j, k;
583
584 for (i = j = k = 0; i < a->isize; ++i)
585 {
586 while (j < b->isize && b->index[j] < a->index[i])
587 {
588 ++j;
589 }
590 if (j == b->isize || b->index[j] != a->index[i])
591 {
592 dest->index[k++] = a->index[i];
593 }
594 }
595 dest->isize = k;
596 }
597
598 /*!
599 * \param[in] a First index group.
600 * \param[in] b Second index group.
601 * \returns Size of the difference \p a - \p b.
602 */
gmx_ana_index_difference_size(gmx_ana_index_t * a,gmx_ana_index_t * b)603 int gmx_ana_index_difference_size(gmx_ana_index_t* a, gmx_ana_index_t* b)
604 {
605 int i, j, k;
606
607 for (i = j = k = 0; i < a->isize; ++i)
608 {
609 while (j < b->isize && b->index[j] < a->index[i])
610 {
611 ++j;
612 }
613 if (j == b->isize || b->index[j] != a->index[i])
614 {
615 ++k;
616 }
617 }
618 return k;
619 }
620
621 /*!
622 * \param[out] dest1 Output group 1 (will equal \p g).
623 * \param[out] dest2 Output group 2 (will equal \p src - \p g).
624 * \param[in] src Group to be partitioned.
625 * \param[in] g One partition.
626 *
627 * \pre \p g is a subset of \p src and both sets are sorted
628 * \pre \p dest1 has allocated storage to store \p src
629 * \post \p dest1 == \p g
630 * \post \p dest2 == \p src - \p g
631 *
632 * No storage should be allocated for \p dest2; after the call,
633 * \p dest2->index points to the memory allocated for \p dest1
634 * (to a part that is not used by \p dest1).
635 *
636 * The calculation can be performed in-place by setting \p dest1 equal to
637 * \p src.
638 */
gmx_ana_index_partition(gmx_ana_index_t * dest1,gmx_ana_index_t * dest2,gmx_ana_index_t * src,gmx_ana_index_t * g)639 void gmx_ana_index_partition(gmx_ana_index_t* dest1, gmx_ana_index_t* dest2, gmx_ana_index_t* src, gmx_ana_index_t* g)
640 {
641 int i, j, k;
642
643 dest2->index = dest1->index + g->isize;
644 dest2->isize = src->isize - g->isize;
645 for (i = g->isize - 1, j = src->isize - 1, k = dest2->isize - 1; i >= 0; --i, --j)
646 {
647 while (j >= 0 && src->index[j] != g->index[i])
648 {
649 dest2->index[k--] = src->index[j--];
650 }
651 }
652 while (j >= 0)
653 {
654 dest2->index[k--] = src->index[j--];
655 }
656 gmx_ana_index_copy(dest1, g, false);
657 }
658
659 /*!
660 * \param[out] dest Output index group (the union of \p a and \p b).
661 * \param[in] a First index group.
662 * \param[in] b Second index group.
663 *
664 * \p a and \p b can have common items.
665 * \p dest can equal \p a or \p b.
666 *
667 * \see gmx_ana_index_merge()
668 */
gmx_ana_index_union(gmx_ana_index_t * dest,gmx_ana_index_t * a,gmx_ana_index_t * b)669 void gmx_ana_index_union(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b)
670 {
671 int dsize;
672 int i, j, k;
673
674 dsize = gmx_ana_index_difference_size(b, a);
675 i = a->isize - 1;
676 j = b->isize - 1;
677 dest->isize = a->isize + dsize;
678 for (k = dest->isize - 1; k >= 0; k--)
679 {
680 if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
681 {
682 dest->index[k] = b->index[j--];
683 }
684 else
685 {
686 if (j >= 0 && a->index[i] == b->index[j])
687 {
688 --j;
689 }
690 dest->index[k] = a->index[i--];
691 }
692 }
693 }
694
gmx_ana_index_union_unsorted(gmx_ana_index_t * dest,gmx_ana_index_t * a,gmx_ana_index_t * b)695 void gmx_ana_index_union_unsorted(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b)
696 {
697 if (gmx_ana_index_check_sorted(b))
698 {
699 gmx_ana_index_union(dest, a, b);
700 }
701 else
702 {
703 gmx_ana_index_t tmp;
704 gmx_ana_index_copy(&tmp, b, true);
705 gmx_ana_index_sort(&tmp);
706 gmx_ana_index_remove_duplicates(&tmp);
707 gmx_ana_index_union(dest, a, &tmp);
708 gmx_ana_index_deinit(&tmp);
709 }
710 }
711
712 /*!
713 * \param[out] dest Output index group (the union of \p a and \p b).
714 * \param[in] a First index group.
715 * \param[in] b Second index group.
716 *
717 * \p a and \p b should not have common items.
718 * \p dest can equal \p a or \p b.
719 *
720 * \see gmx_ana_index_union()
721 */
gmx_ana_index_merge(gmx_ana_index_t * dest,gmx_ana_index_t * a,gmx_ana_index_t * b)722 void gmx_ana_index_merge(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b)
723 {
724 int i, j, k;
725
726 i = a->isize - 1;
727 j = b->isize - 1;
728 dest->isize = a->isize + b->isize;
729 for (k = dest->isize - 1; k >= 0; k--)
730 {
731 if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
732 {
733 dest->index[k] = b->index[j--];
734 }
735 else
736 {
737 dest->index[k] = a->index[i--];
738 }
739 }
740 }
741
742 /********************************************************************
743 * gmx_ana_indexmap_t and related things
744 ********************************************************************/
745
746 /*! \brief
747 * Helper for splitting a sequence of atom indices into groups.
748 *
749 * \param[in] atomIndex Index of the next atom in the sequence.
750 * \param[in] top Topology structure.
751 * \param[in] type Type of group to split into.
752 * \param[in,out] id Variable to receive the next group id.
753 * \returns `true` if \p atomIndex starts a new group in the sequence, i.e.,
754 * if \p *id was changed.
755 *
756 * \p *id should be initialized to `-1` before first call of this function, and
757 * then each atom index in the sequence passed to the function in turn.
758 *
759 * \ingroup module_selection
760 */
next_group_index(int atomIndex,const gmx_mtop_t * top,e_index_t type,int * id)761 static bool next_group_index(int atomIndex, const gmx_mtop_t* top, e_index_t type, int* id)
762 {
763 int prev = *id;
764 switch (type)
765 {
766 case INDEX_ATOM: *id = atomIndex; break;
767 case INDEX_RES:
768 {
769 int resind, molb = 0;
770 mtopGetAtomAndResidueName(top, atomIndex, &molb, nullptr, nullptr, nullptr, &resind);
771 *id = resind;
772 break;
773 }
774 case INDEX_MOL:
775 {
776 int molb = 0;
777 *id = mtopGetMoleculeIndex(top, atomIndex, &molb);
778 break;
779 }
780 case INDEX_UNKNOWN:
781 case INDEX_ALL: *id = 0; break;
782 }
783 return prev != *id;
784 }
785
786 /*!
787 * \param[in,out] t Output block.
788 * \param[in] top Topology structure
789 * (only used if \p type is \ref INDEX_RES or \ref INDEX_MOL, can be NULL
790 * otherwise).
791 * \param[in] g Index group
792 * (can be NULL if \p type is \ref INDEX_UNKNOWN).
793 * \param[in] type Type of partitioning to make.
794 * \param[in] bComplete
795 * If true, the index group is expanded to include any residue/molecule
796 * (depending on \p type) that is partially contained in the group.
797 * If \p type is not INDEX_RES or INDEX_MOL, this has no effect.
798 *
799 * \p m should have been initialized somehow (calloc() is enough).
800 * \p g should be sorted.
801 */
gmx_ana_index_make_block(t_blocka * t,const gmx_mtop_t * top,gmx_ana_index_t * g,e_index_t type,bool bComplete)802 void gmx_ana_index_make_block(t_blocka* t, const gmx_mtop_t* top, gmx_ana_index_t* g, e_index_t type, bool bComplete)
803 {
804 if (type == INDEX_UNKNOWN)
805 {
806 sfree(t->a);
807 srenew(t->index, 2);
808 t->nr = 1;
809 t->nalloc_index = 2;
810 t->index[0] = 0;
811 t->index[1] = 0;
812 t->nra = 0;
813 t->a = nullptr;
814 t->nalloc_a = 0;
815 return;
816 }
817
818 // TODO: Check callers and either check these there as well, or turn these
819 // into exceptions.
820 GMX_RELEASE_ASSERT(top != nullptr || (type != INDEX_RES && type != INDEX_MOL),
821 "Topology must be provided for residue or molecule blocks");
822 GMX_RELEASE_ASSERT(type != INDEX_MOL || top->haveMoleculeIndices,
823 "Molecule information must be present for molecule blocks");
824
825 /* bComplete only does something for INDEX_RES or INDEX_MOL, so turn it
826 * off otherwise. */
827 if (type != INDEX_RES && type != INDEX_MOL)
828 {
829 bComplete = false;
830 }
831 /* Allocate memory for the atom array and fill it unless we are using
832 * completion. */
833 if (bComplete)
834 {
835 t->nra = 0;
836 /* We may allocate some extra memory here because we don't know in
837 * advance how much will be needed. */
838 if (t->nalloc_a < top->natoms)
839 {
840 srenew(t->a, top->natoms);
841 t->nalloc_a = top->natoms;
842 }
843 }
844 else
845 {
846 t->nra = g->isize;
847 if (t->nalloc_a < g->isize)
848 {
849 srenew(t->a, g->isize);
850 t->nalloc_a = g->isize;
851 }
852 if (t->nra > 0)
853 {
854 std::memcpy(t->a, g->index, g->isize * sizeof(*(t->a)));
855 }
856 }
857
858 /* Allocate memory for the block index. We don't know in advance
859 * how much will be needed, so we allocate some extra and free it in the
860 * end. */
861 if (t->nalloc_index < g->isize + 1)
862 {
863 srenew(t->index, g->isize + 1);
864 t->nalloc_index = g->isize + 1;
865 }
866 /* Clear counters */
867 t->nr = 0;
868 int id = -1;
869 int molb = 0;
870 for (int i = 0; i < g->isize; ++i)
871 {
872 const int ai = g->index[i];
873 /* Find the ID number of the atom/residue/molecule corresponding to
874 * the atom. */
875 if (next_group_index(ai, top, type, &id))
876 {
877 /* If this is the first atom in a new block, initialize the block. */
878 if (bComplete)
879 {
880 /* For completion, we first set the start of the block. */
881 t->index[t->nr++] = t->nra;
882 /* And then we find all the atoms that should be included. */
883 switch (type)
884 {
885 case INDEX_RES:
886 {
887 int molnr, atnr_mol;
888 mtopGetMolblockIndex(top, ai, &molb, &molnr, &atnr_mol);
889 const t_atoms& mol_atoms = top->moltype[top->molblock[molb].type].atoms;
890 int last_atom = atnr_mol + 1;
891 const int currentResid = mol_atoms.atom[atnr_mol].resind;
892 while (last_atom < mol_atoms.nr && mol_atoms.atom[last_atom].resind == currentResid)
893 {
894 ++last_atom;
895 }
896 int first_atom = atnr_mol - 1;
897 while (first_atom >= 0 && mol_atoms.atom[first_atom].resind == currentResid)
898 {
899 --first_atom;
900 }
901 const MoleculeBlockIndices& molBlock = top->moleculeBlockIndices[molb];
902 int first_mol_atom = molBlock.globalAtomStart;
903 first_mol_atom += molnr * molBlock.numAtomsPerMolecule;
904 first_atom = first_mol_atom + first_atom + 1;
905 last_atom = first_mol_atom + last_atom - 1;
906 for (int j = first_atom; j <= last_atom; ++j)
907 {
908 t->a[t->nra++] = j;
909 }
910 break;
911 }
912 case INDEX_MOL:
913 {
914 int molnr, atnr_mol;
915 mtopGetMolblockIndex(top, ai, &molb, &molnr, &atnr_mol);
916 const MoleculeBlockIndices& blockIndices = top->moleculeBlockIndices[molb];
917 const int atomStart = blockIndices.globalAtomStart
918 + (id - blockIndices.moleculeIndexStart)
919 * blockIndices.numAtomsPerMolecule;
920 for (int j = 0; j < blockIndices.numAtomsPerMolecule; ++j)
921 {
922 t->a[t->nra++] = atomStart + j;
923 }
924 break;
925 }
926 default: /* Should not be reached */
927 GMX_RELEASE_ASSERT(false, "Unreachable code was reached");
928 break;
929 }
930 }
931 else
932 {
933 /* If not using completion, simply store the start of the block. */
934 t->index[t->nr++] = i;
935 if (type == INDEX_ATOM)
936 {
937 id = -1;
938 }
939 }
940 }
941 }
942 /* Set the end of the last block */
943 t->index[t->nr] = t->nra;
944 /* Free any unnecessary memory */
945 srenew(t->index, t->nr + 1);
946 t->nalloc_index = t->nr + 1;
947 if (bComplete)
948 {
949 srenew(t->a, t->nra);
950 t->nalloc_a = t->nra;
951 }
952 }
953
954 /*!
955 * \param[in] g Index group to check.
956 * \param[in] b Block data to check against.
957 * \returns true if \p g consists of one or more complete blocks from \p b,
958 * false otherwise.
959 *
960 * The atoms in \p g are assumed to be sorted.
961 */
gmx_ana_index_has_full_blocks(const gmx_ana_index_t * g,const gmx::RangePartitioning * b)962 bool gmx_ana_index_has_full_blocks(const gmx_ana_index_t* g, const gmx::RangePartitioning* b)
963 {
964 int i, j, bi;
965
966 i = bi = 0;
967 /* Each round in the loop matches one block */
968 while (i < g->isize)
969 {
970 /* Find the block that begins with the first unmatched atom */
971 while (bi < b->numBlocks() && *b->block(bi).begin() != g->index[i])
972 {
973 ++bi;
974 }
975 /* If not found, or if too large, return */
976 if (bi == b->numBlocks() || i + b->block(bi).size() > g->isize)
977 {
978 return false;
979 }
980 /* Check that the block matches the index */
981 for (j = *b->block(bi).begin(); j < *b->block(bi).end(); ++j, ++i)
982 {
983 if (g->index[i] != j)
984 {
985 return false;
986 }
987 }
988 /* Move the search to the next block */
989 ++bi;
990 }
991 return true;
992 }
993
994 /*!
995 * \param[in] g Index group to check.
996 * \param[in] b Block data to check against.
997 * \returns true if \p g consists of one or more complete blocks from \p b,
998 * false otherwise.
999 *
1000 * The atoms in \p g and \p b->a are assumed to be in the same order.
1001 */
gmx_ana_index_has_full_ablocks(gmx_ana_index_t * g,t_blocka * b)1002 bool gmx_ana_index_has_full_ablocks(gmx_ana_index_t* g, t_blocka* b)
1003 {
1004 int i, j, bi;
1005
1006 i = bi = 0;
1007 /* Each round in the loop matches one block */
1008 while (i < g->isize)
1009 {
1010 /* Find the block that begins with the first unmatched atom */
1011 while (bi < b->nr && b->a[b->index[bi]] != g->index[i])
1012 {
1013 ++bi;
1014 }
1015 /* If not found, or if too large, return */
1016 if (bi == b->nr || i + b->index[bi + 1] - b->index[bi] > g->isize)
1017 {
1018 return false;
1019 }
1020 /* Check that the block matches the index */
1021 for (j = b->index[bi]; j < b->index[bi + 1]; ++j, ++i)
1022 {
1023 if (b->a[j] != g->index[i])
1024 {
1025 return false;
1026 }
1027 }
1028 /* Move the search to the next block */
1029 ++bi;
1030 }
1031 return true;
1032 }
1033
1034 /*!
1035 * \brief Returns if an atom is at a residue boundary.
1036 *
1037 * \param[in] top Topology data.
1038 * \param[in] a Atom index to check, should be -1 <= \p a < top->natoms.
1039 * \param[in,out] molb The molecule block of atom a
1040 * \returns true if atoms \p a and \p a + 1 are in different residues, false otherwise.
1041 */
is_at_residue_boundary(const gmx_mtop_t * top,int a,int * molb)1042 static bool is_at_residue_boundary(const gmx_mtop_t* top, int a, int* molb)
1043 {
1044 if (a == -1 || a + 1 == top->natoms)
1045 {
1046 return true;
1047 }
1048 int resindA;
1049 mtopGetAtomAndResidueName(top, a, molb, nullptr, nullptr, nullptr, &resindA);
1050 int resindAPlusOne;
1051 mtopGetAtomAndResidueName(top, a + 1, molb, nullptr, nullptr, nullptr, &resindAPlusOne);
1052 return resindAPlusOne != resindA;
1053 }
1054
1055 /*!
1056 * \param[in] g Index group to check.
1057 * \param[in] type Block data to check against.
1058 * \param[in] top Topology data.
1059 * \returns true if \p g consists of one or more complete elements of type
1060 * \p type, false otherwise.
1061 *
1062 * \p g is assumed to be sorted, otherwise may return false negatives.
1063 *
1064 * If \p type is \ref INDEX_ATOM, the return value is always true.
1065 * If \p type is \ref INDEX_UNKNOWN or \ref INDEX_ALL, the return value is
1066 * always false.
1067 */
gmx_ana_index_has_complete_elems(gmx_ana_index_t * g,e_index_t type,const gmx_mtop_t * top)1068 bool gmx_ana_index_has_complete_elems(gmx_ana_index_t* g, e_index_t type, const gmx_mtop_t* top)
1069 {
1070 if (g->isize == 0)
1071 {
1072 return true;
1073 }
1074
1075 // TODO: Consider whether unsorted groups need to be supported better.
1076 switch (type)
1077 {
1078 case INDEX_UNKNOWN:
1079 case INDEX_ALL: return false;
1080
1081 case INDEX_ATOM: return true;
1082
1083 case INDEX_RES:
1084 {
1085 int molb = 0;
1086 int aPrev = -1;
1087 for (int i = 0; i < g->isize; ++i)
1088 {
1089 const int a = g->index[i];
1090 // Check if a is consecutive or on a residue boundary
1091 if (a != aPrev + 1)
1092 {
1093 if (!is_at_residue_boundary(top, aPrev, &molb))
1094 {
1095 return false;
1096 }
1097 if (!is_at_residue_boundary(top, a - 1, &molb))
1098 {
1099 return false;
1100 }
1101 }
1102 aPrev = a;
1103 }
1104 GMX_ASSERT(g->isize > 0, "We return above when isize=0");
1105 const int a = g->index[g->isize - 1];
1106 if (!is_at_residue_boundary(top, a, &molb))
1107 {
1108 return false;
1109 }
1110 break;
1111 }
1112
1113 case INDEX_MOL:
1114 {
1115 auto molecules = gmx_mtop_molecules(*top);
1116 return gmx_ana_index_has_full_blocks(g, &molecules);
1117 }
1118 }
1119 return true;
1120 }
1121
1122 /*!
1123 * \param[out] m Output structure.
1124 *
1125 * Any contents of \p m are discarded without freeing.
1126 */
gmx_ana_indexmap_clear(gmx_ana_indexmap_t * m)1127 void gmx_ana_indexmap_clear(gmx_ana_indexmap_t* m)
1128 {
1129 m->type = INDEX_UNKNOWN;
1130 m->refid = nullptr;
1131 m->mapid = nullptr;
1132 m->mapb.nr = 0;
1133 m->mapb.index = nullptr;
1134 m->mapb.nalloc_index = 0;
1135 m->mapb.nra = 0;
1136 m->mapb.a = nullptr;
1137 m->mapb.nalloc_a = 0;
1138 m->orgid = nullptr;
1139 m->b.nr = 0;
1140 m->b.index = nullptr;
1141 m->b.nra = 0;
1142 m->b.a = nullptr;
1143 m->b.nalloc_index = 0;
1144 m->b.nalloc_a = 0;
1145 m->bStatic = true;
1146 }
1147
1148 /*!
1149 * \param[in,out] m Mapping structure.
1150 * \param[in] nr Maximum number of blocks to reserve space for.
1151 * \param[in] isize Maximum number of atoms to reserve space for.
1152 */
gmx_ana_indexmap_reserve(gmx_ana_indexmap_t * m,int nr,int isize)1153 void gmx_ana_indexmap_reserve(gmx_ana_indexmap_t* m, int nr, int isize)
1154 {
1155 if (m->mapb.nalloc_index < nr + 1)
1156 {
1157 srenew(m->refid, nr);
1158 srenew(m->mapid, nr);
1159 srenew(m->orgid, nr);
1160 srenew(m->mapb.index, nr + 1);
1161 srenew(m->b.index, nr + 1);
1162 m->mapb.nalloc_index = nr + 1;
1163 m->b.nalloc_index = nr + 1;
1164 }
1165 if (m->b.nalloc_a < isize)
1166 {
1167 srenew(m->b.a, isize);
1168 m->b.nalloc_a = isize;
1169 }
1170 }
1171
1172 /*!
1173 * \param[in,out] m Mapping structure to initialize.
1174 * \param[in] g Index group to map
1175 * (can be NULL if \p type is \ref INDEX_UNKNOWN).
1176 * \param[in] top Topology structure
1177 * (can be NULL if \p type is not \ref INDEX_RES or \ref INDEX_MOL).
1178 * \param[in] type Type of mapping to construct.
1179 *
1180 * Initializes a new index group mapping.
1181 * The index group provided to gmx_ana_indexmap_update() should always be a
1182 * subset of the \p g given here.
1183 *
1184 * \p m should have been initialized somehow (calloc() is enough).
1185 */
gmx_ana_indexmap_init(gmx_ana_indexmap_t * m,gmx_ana_index_t * g,const gmx_mtop_t * top,e_index_t type)1186 void gmx_ana_indexmap_init(gmx_ana_indexmap_t* m, gmx_ana_index_t* g, const gmx_mtop_t* top, e_index_t type)
1187 {
1188 m->type = type;
1189 gmx_ana_index_make_block(&m->b, top, g, type, false);
1190 gmx_ana_indexmap_reserve(m, m->b.nr, m->b.nra);
1191 int id = -1;
1192 for (int i = 0; i < m->b.nr; ++i)
1193 {
1194 const int ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1195 next_group_index(ii, top, type, &id);
1196 m->refid[i] = i;
1197 m->mapid[i] = id;
1198 m->orgid[i] = id;
1199 }
1200 m->mapb.nr = m->b.nr;
1201 m->mapb.nra = m->b.nra;
1202 m->mapb.a = m->b.a;
1203 std::memcpy(m->mapb.index, m->b.index, (m->b.nr + 1) * sizeof(*(m->mapb.index)));
1204 m->bStatic = true;
1205 }
1206
gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t * m,const gmx_mtop_t * top,e_index_t type)1207 int gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t* m, const gmx_mtop_t* top, e_index_t type)
1208 {
1209 GMX_RELEASE_ASSERT(m->bStatic,
1210 "Changing original IDs is not supported after starting "
1211 "to use the mapping");
1212 GMX_RELEASE_ASSERT(top != nullptr || (type != INDEX_RES && type != INDEX_MOL),
1213 "Topology must be provided for residue or molecule blocks");
1214 // Check that all atoms in each block belong to the same group.
1215 // This is a separate loop for better error handling (no state is modified
1216 // if there is an error.
1217 if (type == INDEX_RES || type == INDEX_MOL)
1218 {
1219 int id = -1;
1220 for (int i = 0; i < m->b.nr; ++i)
1221 {
1222 const int ii = m->b.a[m->b.index[i]];
1223 if (next_group_index(ii, top, type, &id))
1224 {
1225 for (int j = m->b.index[i] + 1; j < m->b.index[i + 1]; ++j)
1226 {
1227 if (next_group_index(m->b.a[j], top, type, &id))
1228 {
1229 std::string message("Grouping into residues/molecules is ambiguous");
1230 GMX_THROW(gmx::InconsistentInputError(message));
1231 }
1232 }
1233 }
1234 }
1235 }
1236 // Do a second loop, where things are actually set.
1237 int id = -1;
1238 int group = -1;
1239 for (int i = 0; i < m->b.nr; ++i)
1240 {
1241 const int ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1242 if (next_group_index(ii, top, type, &id))
1243 {
1244 ++group;
1245 }
1246 m->mapid[i] = group;
1247 m->orgid[i] = group;
1248 }
1249 // Count also the last group.
1250 ++group;
1251 return group;
1252 }
1253
1254 /*!
1255 * \param[in,out] m Mapping structure to initialize.
1256 * \param[in] b Block information to use for data.
1257 *
1258 * Frees some memory that is not necessary for static index group mappings.
1259 * Internal pointers are set to point to data in \p b; it is the responsibility
1260 * of the caller to ensure that the block information matches the contents of
1261 * the mapping.
1262 * After this function has been called, the index group provided to
1263 * gmx_ana_indexmap_update() should always be the same as \p g given here.
1264 *
1265 * This function breaks modularity of the index group mapping interface in an
1266 * ugly way, but allows reducing memory usage of static selections by a
1267 * significant amount.
1268 */
gmx_ana_indexmap_set_static(gmx_ana_indexmap_t * m,t_blocka * b)1269 void gmx_ana_indexmap_set_static(gmx_ana_indexmap_t* m, t_blocka* b)
1270 {
1271 sfree(m->mapid);
1272 sfree(m->mapb.index);
1273 sfree(m->b.index);
1274 sfree(m->b.a);
1275 m->mapb.nalloc_index = 0;
1276 m->mapb.nalloc_a = 0;
1277 m->b.nalloc_index = 0;
1278 m->b.nalloc_a = 0;
1279 m->mapid = m->orgid;
1280 m->mapb.index = b->index;
1281 m->mapb.a = b->a;
1282 m->b.index = b->index;
1283 m->b.a = b->a;
1284 }
1285
1286 /*!
1287 * \param[in,out] dest Destination data structure.
1288 * \param[in] src Source mapping.
1289 * \param[in] bFirst If true, memory is allocated for \p dest and a full
1290 * copy is made; otherwise, only variable parts are copied, and no memory
1291 * is allocated.
1292 *
1293 * \p dest should have been initialized somehow (calloc() is enough).
1294 */
gmx_ana_indexmap_copy(gmx_ana_indexmap_t * dest,gmx_ana_indexmap_t * src,bool bFirst)1295 void gmx_ana_indexmap_copy(gmx_ana_indexmap_t* dest, gmx_ana_indexmap_t* src, bool bFirst)
1296 {
1297 if (bFirst)
1298 {
1299 gmx_ana_indexmap_reserve(dest, src->b.nr, src->b.nra);
1300 dest->type = src->type;
1301 dest->b.nr = src->b.nr;
1302 dest->b.nra = src->b.nra;
1303 std::memcpy(dest->orgid, src->orgid, dest->b.nr * sizeof(*dest->orgid));
1304 std::memcpy(dest->b.index, src->b.index, (dest->b.nr + 1) * sizeof(*dest->b.index));
1305 if (dest->b.nra > 0)
1306 {
1307 std::memcpy(dest->b.a, src->b.a, dest->b.nra * sizeof(*dest->b.a));
1308 }
1309 }
1310 dest->mapb.nr = src->mapb.nr;
1311 dest->mapb.nra = src->mapb.nra;
1312 if (src->mapb.nalloc_a > 0)
1313 {
1314 if (bFirst)
1315 {
1316 snew(dest->mapb.a, src->mapb.nalloc_a);
1317 dest->mapb.nalloc_a = src->mapb.nalloc_a;
1318 }
1319 std::memcpy(dest->mapb.a, src->mapb.a, dest->mapb.nra * sizeof(*dest->mapb.a));
1320 }
1321 else
1322 {
1323 dest->mapb.a = src->mapb.a;
1324 }
1325 std::memcpy(dest->refid, src->refid, dest->mapb.nr * sizeof(*dest->refid));
1326 std::memcpy(dest->mapid, src->mapid, dest->mapb.nr * sizeof(*dest->mapid));
1327 std::memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr + 1) * sizeof(*dest->mapb.index));
1328 dest->bStatic = src->bStatic;
1329 }
1330
1331 /*! \brief
1332 * Helper function to set the source atoms in an index map.
1333 *
1334 * \param[in,out] m Mapping structure.
1335 * \param[in] isize Number of atoms in the \p index array.
1336 * \param[in] index List of atoms.
1337 */
set_atoms(gmx_ana_indexmap_t * m,int isize,int * index)1338 static void set_atoms(gmx_ana_indexmap_t* m, int isize, int* index)
1339 {
1340 m->mapb.nra = isize;
1341 if (m->mapb.nalloc_a == 0)
1342 {
1343 m->mapb.a = index;
1344 }
1345 else
1346 {
1347 for (int i = 0; i < isize; ++i)
1348 {
1349 m->mapb.a[i] = index[i];
1350 }
1351 }
1352 }
1353
1354 /*!
1355 * \param[in,out] m Mapping structure.
1356 * \param[in] g Current index group.
1357 * \param[in] bMaskOnly true if the unused blocks should be masked with
1358 * -1 instead of removing them.
1359 *
1360 * Updates the index group mapping with the new index group \p g.
1361 *
1362 * \see gmx_ana_indexmap_t
1363 */
gmx_ana_indexmap_update(gmx_ana_indexmap_t * m,gmx_ana_index_t * g,bool bMaskOnly)1364 void gmx_ana_indexmap_update(gmx_ana_indexmap_t* m, gmx_ana_index_t* g, bool bMaskOnly)
1365 {
1366 int i, j, bi, bj;
1367
1368 /* Process the simple cases first */
1369 if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
1370 {
1371 return;
1372 }
1373 if (m->type == INDEX_ALL)
1374 {
1375 set_atoms(m, g->isize, g->index);
1376 if (m->b.nr > 0)
1377 {
1378 m->mapb.index[1] = g->isize;
1379 }
1380 return;
1381 }
1382 /* Reset the reference IDs and mapping if necessary */
1383 const bool bToFull = (g->isize == m->b.nra);
1384 const bool bWasFull = (m->mapb.nra == m->b.nra);
1385 if (bToFull || bMaskOnly)
1386 {
1387 if (!m->bStatic)
1388 {
1389 for (bj = 0; bj < m->b.nr; ++bj)
1390 {
1391 m->refid[bj] = bj;
1392 }
1393 }
1394 if (!bWasFull)
1395 {
1396 for (bj = 0; bj < m->b.nr; ++bj)
1397 {
1398 m->mapid[bj] = m->orgid[bj];
1399 }
1400 for (bj = 0; bj <= m->b.nr; ++bj)
1401 {
1402 m->mapb.index[bj] = m->b.index[bj];
1403 }
1404 }
1405 set_atoms(m, m->b.nra, m->b.a);
1406 m->mapb.nr = m->b.nr;
1407 }
1408 /* Exit immediately if the group is static */
1409 if (bToFull)
1410 {
1411 m->bStatic = true;
1412 return;
1413 }
1414
1415 if (bMaskOnly)
1416 {
1417 for (i = j = bj = 0; i < g->isize; ++i, ++j)
1418 {
1419 /* Find the next atom in the block */
1420 while (m->b.a[j] != g->index[i])
1421 {
1422 ++j;
1423 }
1424 /* Mark blocks that did not contain any atoms */
1425 while (bj < m->b.nr && m->b.index[bj + 1] <= j)
1426 {
1427 m->refid[bj++] = -1;
1428 }
1429 /* Advance the block index if we have reached the next block */
1430 if (m->b.index[bj] <= j)
1431 {
1432 ++bj;
1433 }
1434 }
1435 /* Mark the last blocks as not accessible */
1436 while (bj < m->b.nr)
1437 {
1438 m->refid[bj++] = -1;
1439 }
1440 }
1441 else
1442 {
1443 set_atoms(m, g->isize, g->index);
1444 for (i = j = bi = 0, bj = -1; i < g->isize; ++i)
1445 {
1446 /* Find the next atom in the block */
1447 while (m->b.a[j] != g->index[i])
1448 {
1449 ++j;
1450 }
1451 /* If we have reached a new block, add it */
1452 if (m->b.index[bj + 1] <= j)
1453 {
1454 /* Skip any blocks in between */
1455 while (bj < m->b.nr && m->b.index[bj + 1] <= j)
1456 {
1457 ++bj;
1458 }
1459 m->refid[bi] = bj;
1460 m->mapid[bi] = m->orgid[bj];
1461 m->mapb.index[bi] = i;
1462 bi++;
1463 }
1464 }
1465 /* Update the number of blocks */
1466 m->mapb.index[bi] = g->isize;
1467 m->mapb.nr = bi;
1468 }
1469 m->bStatic = false;
1470 }
1471
1472 /*!
1473 * \param[in,out] m Mapping structure to free.
1474 *
1475 * All the memory allocated for the mapping structure is freed, and
1476 * the pointers set to NULL.
1477 * The pointer \p m is not freed.
1478 */
gmx_ana_indexmap_deinit(gmx_ana_indexmap_t * m)1479 void gmx_ana_indexmap_deinit(gmx_ana_indexmap_t* m)
1480 {
1481 sfree(m->refid);
1482 if (m->mapid != m->orgid)
1483 {
1484 sfree(m->mapid);
1485 }
1486 if (m->mapb.nalloc_index > 0)
1487 {
1488 sfree(m->mapb.index);
1489 }
1490 if (m->mapb.nalloc_a > 0)
1491 {
1492 sfree(m->mapb.a);
1493 }
1494 sfree(m->orgid);
1495 if (m->b.nalloc_index > 0)
1496 {
1497 sfree(m->b.index);
1498 }
1499 if (m->b.nalloc_a > 0)
1500 {
1501 sfree(m->b.a);
1502 }
1503 gmx_ana_indexmap_clear(m);
1504 }
1505