src/gmxlib/trajana/indexutil.c

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