src/gmxlib/trajana/indexutil.c

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