src/gromacs/selection/indexutil.cpp

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