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