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