src/gromacs/selection/indexutil.cpp

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