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