3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
33 * Implements functions in indexutil.h.
35 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
36 * \ingroup module_selection
42 #include "gromacs/legacyheaders/index.h"
43 #include "gromacs/legacyheaders/gmx_fatal.h"
44 #include "gromacs/legacyheaders/smalloc.h"
45 #include "gromacs/legacyheaders/string2.h"
46 #include "gromacs/legacyheaders/typedefs.h"
48 #include "gromacs/selection/indexutil.h"
50 /********************************************************************
51 * gmx_ana_indexgrps_t functions
52 ********************************************************************/
55 * Stores a set of index groups.
57 struct gmx_ana_indexgrps_t
59 /** Number of index groups. */
61 /** Array of index groups. */
66 * \param[out] g Index group structure.
67 * \param[in] ngrps Number of groups for which memory is allocated.
70 gmx_ana_indexgrps_alloc(gmx_ana_indexgrps_t **g, int ngrps)
78 * \param[out] g Index group structure.
79 * \param[in] top Topology structure.
80 * \param[in] fnm File name for the index file.
81 * Memory is automatically allocated.
83 * One or both of \p top or \p fnm can be NULL.
84 * If \p top is NULL, an index file is required and the groups are read
85 * from the file (uses Gromacs routine init_index()).
86 * If \p fnm is NULL, default groups are constructed based on the
87 * topology (uses Gromacs routine analyse()).
88 * If both are null, the index group structure is initialized empty.
91 gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, t_topology *top,
94 t_blocka *block = NULL;
100 block = init_index(fnm, &names);
104 block = new_blocka();
105 analyse(&top->atoms, block, &names, FALSE, FALSE);
115 gmx_ana_indexgrps_alloc(g, block->nr);
116 for (i = 0; i < block->nr; ++i)
118 gmx_ana_index_t *grp = &(*g)->g[i];
120 grp->isize = block->index[i+1] - block->index[i];
121 snew(grp->index, grp->isize);
122 for (j = 0; j < grp->isize; ++j)
124 grp->index[j] = block->a[block->index[i]+j];
126 grp->name = names[i];
127 grp->nalloc_index = grp->isize;
136 * \param[in] g Index groups structure.
138 * The pointer \p g is invalid after the call.
141 gmx_ana_indexgrps_free(gmx_ana_indexgrps_t *g)
150 for (i = 0; i < g->nr; ++i)
152 gmx_ana_index_deinit(&g->g[i]);
161 * \param[out] dest Destination index groups.
162 * \param[in] src Source index groups.
164 * A deep copy is made for all fields, including the group names.
167 gmx_ana_indexgrps_clone(gmx_ana_indexgrps_t **dest, gmx_ana_indexgrps_t *src)
171 gmx_ana_indexgrps_alloc(dest, src->nr);
172 for (g = 0; g < src->nr; ++g)
174 gmx_ana_index_copy(&(*dest)->g[g], &src->g[g], true);
179 * \param[out] g Index group structure.
180 * \returns true if \p g is empty, i.e., has 0 index groups.
183 gmx_ana_indexgrps_is_empty(gmx_ana_indexgrps_t *g)
189 * \param[in] g Index groups structure.
190 * \param[in] n Index group number to get.
191 * \returns Pointer to the \p n'th index group in \p g.
193 * The returned pointer should not be freed.
196 gmx_ana_indexgrps_get_grp(gmx_ana_indexgrps_t *g, int n)
198 if (n < 0 || n >= g->nr)
206 * \param[out] dest Output structure.
207 * \param[in] src Input index groups.
208 * \param[in] n Number of the group to extract.
209 * \returns true if \p n is a valid group in \p src, false otherwise.
212 gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, int n)
214 if (n < 0 || n >= src->nr)
220 gmx_ana_index_copy(dest, &src->g[n], true);
225 * \param[out] dest Output structure.
226 * \param[in] src Input index groups.
227 * \param[in] name Name (or part of the name) of the group to extract.
228 * \returns true if \p name is a valid group in \p src, false otherwise.
230 * Uses the Gromacs routine find_group() to find the actual group;
231 * the comparison is case-insensitive.
234 gmx_ana_indexgrps_find(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, char *name)
239 snew(names, src->nr);
240 for (i = 0; i < src->nr; ++i)
242 names[i] = src->g[i].name;
244 i = find_group(name, src->nr, names);
252 return gmx_ana_indexgrps_extract(dest, src, i);
256 * \param[in] fp Where to print the output.
257 * \param[in] g Index groups to print.
258 * \param[in] maxn Maximum number of indices to print
259 * (-1 = print all, 0 = print only names).
262 gmx_ana_indexgrps_print(FILE *fp, gmx_ana_indexgrps_t *g, int maxn)
266 for (i = 0; i < g->nr; ++i)
268 fprintf(fp, " %2d: ", i);
269 gmx_ana_index_dump(fp, &g->g[i], i, maxn);
273 /********************************************************************
274 * gmx_ana_index_t functions
275 ********************************************************************/
278 * \param[in,out] g Index group structure.
279 * \param[in] isize Maximum number of atoms to reserve space for.
282 gmx_ana_index_reserve(gmx_ana_index_t *g, int isize)
284 if (g->nalloc_index < isize)
286 srenew(g->index, isize);
287 g->nalloc_index = isize;
292 * \param[in,out] g Index group structure.
294 * Resizes the memory allocated for holding the indices such that the
295 * current contents fit.
298 gmx_ana_index_squeeze(gmx_ana_index_t *g)
300 srenew(g->index, g->isize);
301 g->nalloc_index = g->isize;
305 * \param[out] g Output structure.
307 * Any contents of \p g are discarded without freeing.
310 gmx_ana_index_clear(gmx_ana_index_t *g)
319 * \param[out] g Output structure.
320 * \param[in] isize Number of atoms in the new group.
321 * \param[in] index Array of \p isize atoms (can be NULL if \p isize is 0).
322 * \param[in] name Name for the new group (can be NULL).
323 * \param[in] nalloc Number of elements allocated for \p index
324 * (if 0, \p index is not freed in gmx_ana_index_deinit())
326 * No copy if \p index is made.
329 gmx_ana_index_set(gmx_ana_index_t *g, int isize, atom_id *index, char *name,
335 g->nalloc_index = nalloc;
339 * \param[out] g Output structure.
340 * \param[in] natoms Number of atoms.
341 * \param[in] name Name for the new group (can be NULL).
344 gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms, char *name)
349 snew(g->index, natoms);
350 for (i = 0; i < natoms; ++i)
355 g->nalloc_index = natoms;
359 * \param[in] g Index group structure.
361 * The pointer \p g is not freed.
364 gmx_ana_index_deinit(gmx_ana_index_t *g)
366 if (g->nalloc_index > 0)
371 gmx_ana_index_clear(g);
375 * \param[out] dest Destination index group.
376 * \param[in] src Source index group.
377 * \param[in] bAlloc If true, memory is allocated at \p dest; otherwise,
378 * it is assumed that enough memory has been allocated for index.
380 * A deep copy of the name is only made if \p bAlloc is true.
383 gmx_ana_index_copy(gmx_ana_index_t *dest, gmx_ana_index_t *src, bool bAlloc)
385 dest->isize = src->isize;
390 snew(dest->index, dest->isize);
391 dest->nalloc_index = dest->isize;
393 memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
395 if (bAlloc && src->name)
397 dest->name = strdup(src->name);
399 else if (bAlloc || src->name)
401 dest->name = src->name;
406 * \param[in] fp Where to print the output.
407 * \param[in] g Index group to print.
408 * \param[in] i Group number to use if the name is NULL.
409 * \param[in] maxn Maximum number of indices to print (-1 = print all).
412 gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int i, int maxn)
418 fprintf(fp, "\"%s\"", g->name);
422 fprintf(fp, "Group %d", i+1);
424 fprintf(fp, " (%d atoms)", g->isize);
429 if (maxn >= 0 && n > maxn)
433 for (j = 0; j < n; ++j)
435 fprintf(fp, " %d", g->index[j]+1);
446 * \param[in] g Input index group.
447 * \param[in] natoms Number of atoms to check against.
449 * If any atom index in the index group is less than zero or >= \p natoms,
450 * gmx_fatal() is called.
453 gmx_ana_index_check(gmx_ana_index_t *g, int natoms)
457 for (j = 0; j < g->isize; ++j)
459 if (g->index[j] >= natoms)
461 gmx_fatal(FARGS,"Atom index (%d) in index group %s (%d atoms) "
462 "larger than number of atoms in trajectory (%d atoms)",
463 g->index[j], g->name, g->isize, natoms);
465 else if (g->index[j] < 0)
467 gmx_fatal(FARGS,"Atom index (%d) in index group %s (%d atoms) "
469 g->index[j], g->name, g->isize);
475 * \param[in] g Index group to check.
476 * \returns true if the index group is sorted and has no duplicates,
480 gmx_ana_index_check_sorted(gmx_ana_index_t *g)
484 for (i = 0; i < g->isize-1; ++i)
486 if (g->index[i+1] <= g->index[i])
494 /********************************************************************
496 ********************************************************************/
498 /** Helper function for gmx_ana_index_sort(). */
500 cmp_atomid(const void *a, const void *b)
502 if (*(atom_id *)a < *(atom_id *)b) return -1;
503 if (*(atom_id *)a > *(atom_id *)b) return 1;
508 * \param[in,out] g Index group to be sorted.
511 gmx_ana_index_sort(gmx_ana_index_t *g)
513 qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
517 * \param[in] a Index group to check.
518 * \param[in] b Index group to check.
519 * \returns true if \p a and \p b are equal, false otherwise.
522 gmx_ana_index_equals(gmx_ana_index_t *a, gmx_ana_index_t *b)
526 if (a->isize != b->isize)
530 for (i = 0; i < a->isize; ++i)
532 if (a->index[i] != b->index[i])
541 * \param[in] a Index group to check against.
542 * \param[in] b Index group to check.
543 * \returns true if \p b is contained in \p a,
546 * If the elements are not in the same order in both groups, the function
547 * fails. However, the groups do not need to be sorted.
550 gmx_ana_index_contains(gmx_ana_index_t *a, gmx_ana_index_t *b)
554 for (i = j = 0; j < b->isize; ++i, ++j) {
555 while (i < a->isize && a->index[i] != b->index[j])
568 * \param[out] dest Output index group (the intersection of \p a and \p b).
569 * \param[in] a First index group.
570 * \param[in] b Second index group.
572 * \p dest can be the same as \p a or \p b.
575 gmx_ana_index_intersection(gmx_ana_index_t *dest,
576 gmx_ana_index_t *a, gmx_ana_index_t *b)
580 for (i = j = k = 0; i < a->isize && j < b->isize; ++i) {
581 while (j < b->isize && b->index[j] < a->index[i])
585 if (j < b->isize && b->index[j] == a->index[i])
587 dest->index[k++] = b->index[j++];
594 * \param[out] dest Output index group (the difference \p a - \p b).
595 * \param[in] a First index group.
596 * \param[in] b Second index group.
598 * \p dest can equal \p a, but not \p b.
601 gmx_ana_index_difference(gmx_ana_index_t *dest,
602 gmx_ana_index_t *a, gmx_ana_index_t *b)
606 for (i = j = k = 0; i < a->isize; ++i)
608 while (j < b->isize && b->index[j] < a->index[i])
612 if (j == b->isize || b->index[j] != a->index[i])
614 dest->index[k++] = a->index[i];
621 * \param[in] a First index group.
622 * \param[in] b Second index group.
623 * \returns Size of the difference \p a - \p b.
626 gmx_ana_index_difference_size(gmx_ana_index_t *a, gmx_ana_index_t *b)
630 for (i = j = k = 0; i < a->isize; ++i)
632 while (j < b->isize && b->index[j] < a->index[i])
636 if (j == b->isize || b->index[j] != a->index[i])
645 * \param[out] dest1 Output group 1 (will equal \p g).
646 * \param[out] dest2 Output group 2 (will equal \p src - \p g).
647 * \param[in] src Group to be partitioned.
648 * \param[in] g One partition.
650 * \pre \p g is a subset of \p src and both sets are sorted
651 * \pre \p dest1 has allocated storage to store \p src
652 * \post \p dest1 == \p g
653 * \post \p dest2 == \p src - \p g
655 * No storage should be allocated for \p dest2; after the call,
656 * \p dest2->index points to the memory allocated for \p dest1
657 * (to a part that is not used by \p dest1).
659 * The calculation can be performed in-place by setting \p dest1 equal to
663 gmx_ana_index_partition(gmx_ana_index_t *dest1, gmx_ana_index_t *dest2,
664 gmx_ana_index_t *src, gmx_ana_index_t *g)
668 dest2->index = dest1->index + g->isize;
669 dest2->isize = src->isize - g->isize;
670 for (i = g->isize-1, j = src->isize-1, k = dest2->isize-1; i >= 0; --i, --j)
672 while (j >= 0 && src->index[j] != g->index[i])
674 dest2->index[k--] = src->index[j--];
679 dest2->index[k--] = src->index[j--];
681 gmx_ana_index_copy(dest1, g, false);
685 * \param[out] dest Output index group (the union of \p a and \p b).
686 * \param[in] a First index group.
687 * \param[in] b Second index group.
689 * \p a and \p b can have common items.
690 * \p dest can equal \p a or \p b.
692 * \see gmx_ana_index_merge()
695 gmx_ana_index_union(gmx_ana_index_t *dest,
696 gmx_ana_index_t *a, gmx_ana_index_t *b)
701 dsize = gmx_ana_index_difference_size(b, a);
704 dest->isize = a->isize + dsize;
705 for (k = dest->isize - 1; k >= 0; k--)
707 if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
709 dest->index[k] = b->index[j--];
713 if (j >= 0 && a->index[i] == b->index[j])
717 dest->index[k] = a->index[i--];
723 * \param[out] dest Output index group (the union of \p a and \p b).
724 * \param[in] a First index group.
725 * \param[in] b Second index group.
727 * \p a and \p b should not have common items.
728 * \p dest can equal \p a or \p b.
730 * \see gmx_ana_index_union()
733 gmx_ana_index_merge(gmx_ana_index_t *dest,
734 gmx_ana_index_t *a, gmx_ana_index_t *b)
740 dest->isize = a->isize + b->isize;
741 for (k = dest->isize - 1; k >= 0; k--)
743 if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
745 dest->index[k] = b->index[j--];
749 dest->index[k] = a->index[i--];
754 /********************************************************************
755 * gmx_ana_indexmap_t and related things
756 ********************************************************************/
759 * \param[in,out] t Output block.
760 * \param[in] top Topology structure
761 * (only used if \p type is \ref INDEX_RES or \ref INDEX_MOL, can be NULL
763 * \param[in] g Index group
764 * (can be NULL if \p type is \ref INDEX_UNKNOWN).
765 * \param[in] type Type of partitioning to make.
766 * \param[in] bComplete
767 * If true, the index group is expanded to include any residue/molecule
768 * (depending on \p type) that is partially contained in the group.
769 * If \p type is not INDEX_RES or INDEX_MOL, this has no effect.
771 * \p m should have been initialized somehow (calloc() is enough) unless
772 * \p type is INDEX_UNKNOWN.
773 * \p g should be sorted.
776 gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
777 e_index_t type, bool bComplete)
782 if (type == INDEX_UNKNOWN)
795 /* bComplete only does something for INDEX_RES or INDEX_MOL, so turn it
797 if (type != INDEX_RES && type != INDEX_MOL)
801 /* Allocate memory for the atom array and fill it unless we are using
806 /* We may allocate some extra memory here because we don't know in
807 * advance how much will be needed. */
808 if (t->nalloc_a < top->atoms.nr)
810 srenew(t->a, top->atoms.nr);
811 t->nalloc_a = top->atoms.nr;
817 if (t->nalloc_a < g->isize)
819 srenew(t->a, g->isize);
820 t->nalloc_a = g->isize;
822 memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
825 /* Allocate memory for the block index. We don't know in advance
826 * how much will be needed, so we allocate some extra and free it in the
828 if (t->nalloc_index < g->isize + 1)
830 srenew(t->index, g->isize + 1);
831 t->nalloc_index = g->isize + 1;
835 j = 0; /* j is used by residue completion for the first atom not stored */
837 for (i = 0; i < g->isize; ++i)
840 /* Find the ID number of the atom/residue/molecule corresponding to
848 id = top->atoms.atom[ai].resind;
851 while (ai >= top->mols.index[id+1])
856 case INDEX_UNKNOWN: /* Should not occur */
861 /* If this is the first atom in a new block, initialize the block. */
866 /* For completion, we first set the start of the block. */
867 t->index[t->nr++] = t->nra;
868 /* And then we find all the atoms that should be included. */
872 while (top->atoms.atom[j].resind != id)
876 while (j < top->atoms.nr && top->atoms.atom[j].resind == id)
884 for (j = top->mols.index[id]; j < top->mols.index[id+1]; ++j)
890 default: /* Should not be reached */
891 gmx_bug("internal error");
897 /* If not using completion, simply store the start of the block. */
898 t->index[t->nr++] = i;
903 /* Set the end of the last block */
904 t->index[t->nr] = t->nra;
905 /* Free any unnecessary memory */
906 srenew(t->index, t->nr+1);
907 t->nalloc_index = t->nr+1;
910 srenew(t->a, t->nra);
911 t->nalloc_a = t->nra;
916 * \param[in] g Index group to check.
917 * \param[in] b Block data to check against.
918 * \returns true if \p g consists of one or more complete blocks from \p b,
921 * The atoms in \p g are assumed to be sorted.
924 gmx_ana_index_has_full_blocks(gmx_ana_index_t *g, t_block *b)
929 /* Each round in the loop matches one block */
932 /* Find the block that begins with the first unmatched atom */
933 while (bi < b->nr && b->index[bi] != g->index[i])
937 /* If not found, or if too large, return */
938 if (bi == b->nr || i + b->index[bi+1] - b->index[bi] > g->isize)
942 /* Check that the block matches the index */
943 for (j = b->index[bi]; j < b->index[bi+1]; ++j, ++i)
945 if (g->index[i] != j)
950 /* Move the search to the next block */
957 * \param[in] g Index group to check.
958 * \param[in] b Block data to check against.
959 * \returns true if \p g consists of one or more complete blocks from \p b,
962 * The atoms in \p g and \p b->a are assumed to be in the same order.
965 gmx_ana_index_has_full_ablocks(gmx_ana_index_t *g, t_blocka *b)
970 /* Each round in the loop matches one block */
973 /* Find the block that begins with the first unmatched atom */
974 while (bi < b->nr && b->a[b->index[bi]] != g->index[i])
978 /* If not found, or if too large, return */
979 if (bi == b->nr || i + b->index[bi+1] - b->index[bi] > g->isize)
983 /* Check that the block matches the index */
984 for (j = b->index[bi]; j < b->index[bi+1]; ++j, ++i)
986 if (b->a[j] != g->index[i])
991 /* Move the search to the next block */
998 * \param[in] g Index group to check.
999 * \param[in] type Block data to check against.
1000 * \param[in] top Topology data.
1001 * \returns true if \p g consists of one or more complete elements of type
1002 * \p type, false otherwise.
1004 * If \p type is \ref INDEX_ATOM, the return value is always true.
1005 * If \p type is \ref INDEX_UNKNOWN or \ref INDEX_ALL, the return value is
1009 gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
1027 for (i = 0; i < g->isize; ++i)
1030 id = top->atoms.atom[ai].resind;
1033 if (ai > 0 && top->atoms.atom[ai-1].resind == id)
1037 if (i > 0 && g->index[i-1] < top->atoms.nr - 1
1038 && top->atoms.atom[g->index[i-1]+1].resind == prev)
1045 if (g->index[i-1] < top->atoms.nr - 1
1046 && top->atoms.atom[g->index[i-1]+1].resind == prev)
1054 return gmx_ana_index_has_full_blocks(g, &top->mols);
1060 * \param[out] m Output structure.
1062 * Any contents of \p m are discarded without freeing.
1065 gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
1067 m->type = INDEX_UNKNOWN;
1072 m->mapb.index = NULL;
1073 m->mapb.nalloc_index = 0;
1079 m->b.nalloc_index = 0;
1082 m->bMapStatic = true;
1086 * \param[in,out] m Mapping structure.
1087 * \param[in] nr Maximum number of blocks to reserve space for.
1088 * \param[in] isize Maximum number of atoms to reserve space for.
1091 gmx_ana_indexmap_reserve(gmx_ana_indexmap_t *m, int nr, int isize)
1093 if (m->mapb.nalloc_index < nr + 1)
1095 srenew(m->refid, nr);
1096 srenew(m->mapid, nr);
1097 srenew(m->orgid, nr);
1098 srenew(m->mapb.index, nr + 1);
1099 srenew(m->b.index, nr + 1);
1100 m->mapb.nalloc_index = nr + 1;
1101 m->b.nalloc_index = nr + 1;
1103 if (m->b.nalloc_a < isize)
1105 srenew(m->b.a, isize);
1106 m->b.nalloc_a = isize;
1111 * \param[in,out] m Mapping structure to initialize.
1112 * \param[in] g Index group to map
1113 * (can be NULL if \p type is \ref INDEX_UNKNOWN).
1114 * \param[in] top Topology structure
1115 * (can be NULL if \p type is not \ref INDEX_RES or \ref INDEX_MOL).
1116 * \param[in] type Type of mapping to construct.
1118 * Initializes a new index group mapping.
1119 * The index group provided to gmx_ana_indexmap_update() should always be a
1120 * subset of the \p g given here.
1122 * \p m should have been initialized somehow (calloc() is enough).
1125 gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
1126 t_topology *top, e_index_t type)
1131 gmx_ana_index_make_block(&m->b, top, g, type, false);
1132 gmx_ana_indexmap_reserve(m, m->b.nr, m->b.nra);
1134 for (i = mi = 0; i < m->nr; ++i)
1136 ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1143 m->orgid[i] = top->atoms.atom[ii].resind;
1146 while (top->mols.index[mi+1] <= ii)
1158 for (i = 0; i < m->nr; ++i)
1161 m->mapid[i] = m->orgid[i];
1164 memcpy(m->mapb.index, m->b.index, (m->nr+1)*sizeof(*(m->mapb.index)));
1166 m->bMapStatic = true;
1170 * \param[in,out] m Mapping structure to initialize.
1171 * \param[in] b Block information to use for data.
1173 * Frees some memory that is not necessary for static index group mappings.
1174 * Internal pointers are set to point to data in \p b; it is the responsibility
1175 * of the caller to ensure that the block information matches the contents of
1177 * After this function has been called, the index group provided to
1178 * gmx_ana_indexmap_update() should always be the same as \p g given here.
1180 * This function breaks modularity of the index group mapping interface in an
1181 * ugly way, but allows reducing memory usage of static selections by a
1182 * significant amount.
1185 gmx_ana_indexmap_set_static(gmx_ana_indexmap_t *m, t_blocka *b)
1188 m->mapid = m->orgid;
1190 m->b.nalloc_index = 0;
1191 m->b.index = b->index;
1192 sfree(m->mapb.index);
1193 m->mapb.nalloc_index = 0;
1194 m->mapb.index = m->b.index;
1201 * \param[in,out] dest Destination data structure.
1202 * \param[in] src Source mapping.
1203 * \param[in] bFirst If true, memory is allocated for \p dest and a full
1204 * copy is made; otherwise, only variable parts are copied, and no memory
1207 * \p dest should have been initialized somehow (calloc() is enough).
1210 gmx_ana_indexmap_copy(gmx_ana_indexmap_t *dest, gmx_ana_indexmap_t *src, bool bFirst)
1214 gmx_ana_indexmap_reserve(dest, src->b.nr, src->b.nra);
1215 dest->type = src->type;
1216 dest->b.nr = src->b.nr;
1217 dest->b.nra = src->b.nra;
1218 memcpy(dest->orgid, src->orgid, dest->b.nr*sizeof(*dest->orgid));
1219 memcpy(dest->b.index, src->b.index, (dest->b.nr+1)*sizeof(*dest->b.index));
1220 memcpy(dest->b.a, src->b.a, dest->b.nra*sizeof(*dest->b.a));
1223 dest->mapb.nr = src->mapb.nr;
1224 memcpy(dest->refid, src->refid, dest->nr*sizeof(*dest->refid));
1225 memcpy(dest->mapid, src->mapid, dest->nr*sizeof(*dest->mapid));
1226 memcpy(dest->mapb.index, src->mapb.index,(dest->mapb.nr+1)*sizeof(*dest->mapb.index));
1227 dest->bStatic = src->bStatic;
1228 dest->bMapStatic = src->bMapStatic;
1232 * \param[in,out] m Mapping structure.
1233 * \param[in] g Current index group.
1234 * \param[in] bMaskOnly true if the unused blocks should be masked with
1235 * -1 instead of removing them.
1237 * Updates the index group mapping with the new index group \p g.
1239 * \see gmx_ana_indexmap_t
1242 gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
1248 /* Process the simple cases first */
1249 if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
1253 if (m->type == INDEX_ALL)
1257 m->mapb.index[1] = g->isize;
1261 /* Reset the reference IDs and mapping if necessary */
1262 bStatic = (g->isize == m->b.nra && m->nr == m->b.nr);
1263 if (bStatic || bMaskOnly)
1267 for (bj = 0; bj < m->b.nr; ++bj)
1274 for (bj = 0; bj < m->b.nr; ++bj)
1276 m->mapid[bj] = m->orgid[bj];
1278 for (bj = 0; bj <= m->b.nr; ++bj)
1280 m->mapb.index[bj] = m->b.index[bj];
1282 m->bMapStatic = true;
1285 /* Exit immediately if the group is static */
1295 for (i = j = bj = 0; i < g->isize; ++i, ++j)
1297 /* Find the next atom in the block */
1298 while (m->b.a[j] != g->index[i])
1302 /* Mark blocks that did not contain any atoms */
1303 while (bj < m->b.nr && m->b.index[bj+1] <= j)
1305 m->refid[bj++] = -1;
1307 /* Advance the block index if we have reached the next block */
1308 if (m->b.index[bj] <= j)
1313 /* Mark the last blocks as not accessible */
1314 while (bj < m->b.nr)
1316 m->refid[bj++] = -1;
1321 for (i = j = bi = 0, bj = -1; i < g->isize; ++i)
1323 /* Find the next atom in the block */
1324 while (m->b.a[j] != g->index[i])
1328 /* If we have reached a new block, add it */
1329 if (m->b.index[bj+1] <= j)
1331 /* Skip any blocks in between */
1332 while (bj < m->b.nr && m->b.index[bj+1] <= j)
1337 m->mapid[bi] = m->orgid[bj];
1338 m->mapb.index[bi] = i;
1342 /* Update the number of blocks */
1343 m->mapb.index[bi] = g->isize;
1345 m->bMapStatic = false;
1352 * \param[in,out] m Mapping structure to free.
1354 * All the memory allocated for the mapping structure is freed, and
1355 * the pointers set to NULL.
1356 * The pointer \p m is not freed.
1359 gmx_ana_indexmap_deinit(gmx_ana_indexmap_t *m)
1362 if (m->mapid != m->orgid)
1366 if (m->mapb.nalloc_index > 0)
1368 sfree(m->mapb.index);
1371 if (m->b.nalloc_index > 0)
1375 if (m->b.nalloc_a > 0)
1379 gmx_ana_indexmap_clear(m);