src/gromacs/selection/indexutil.cpp

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