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