src/gromacs/selection/selelem.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014, 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 selelem.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include <cstring>
  43
  44 #include "gromacs/selection/indexutil.h"
  45 #include "gromacs/selection/poscalc.h"
  46 #include "gromacs/selection/position.h"
  47 #include "gromacs/selection/selmethod.h"
  48 #include "gromacs/utility/exceptions.h"
  49 #include "gromacs/utility/gmxassert.h"
  50 #include "gromacs/utility/smalloc.h"
  51 #include "gromacs/utility/stringutil.h"
  52
  53 #include "keywords.h"
  54 #include "mempool.h"
  55 #include "selelem.h"
  56
  57 /*!
  58  * \param[in] sel Selection for which the string is requested
  59  * \returns   Pointer to a string that corresponds to \p sel->type.
  60  *
  61  * The return value points to a string constant and should not be \p free'd.
  62  *
  63  * The function returns NULL if \p sel->type is not one of the valid values.
  64  */
  65 const char *
  66 _gmx_selelem_type_str(const gmx::SelectionTreeElement &sel)
  67 {
  68     switch (sel.type)
  69     {
  70         case SEL_CONST:      return "CONST";
  71         case SEL_EXPRESSION: return "EXPR";
  72         case SEL_BOOLEAN:    return "BOOL";
  73         case SEL_ARITHMETIC: return "ARITH";
  74         case SEL_ROOT:       return "ROOT";
  75         case SEL_SUBEXPR:    return "SUBEXPR";
  76         case SEL_SUBEXPRREF: return "REF";
  77         case SEL_GROUPREF:   return "GROUPREF";
  78         case SEL_MODIFIER:   return "MODIFIER";
  79     }
  80     return NULL;
  81 }
  82
  83 /*!
  84  * \param[in] val Value structore for which the string is requested.
  85  * \returns   Pointer to a string that corresponds to \p val->type,
  86  *   NULL if the type value is invalid.
  87  *
  88  * The return value points to a string constant and should not be \p free'd.
  89  */
  90 const char *
  91 _gmx_sel_value_type_str(const gmx_ana_selvalue_t *val)
  92 {
  93     switch (val->type)
  94     {
  95         case NO_VALUE:       return "NONE";
  96         case INT_VALUE:      return "INT";
  97         case REAL_VALUE:     return "REAL";
  98         case STR_VALUE:      return "STR";
  99         case POS_VALUE:      return "VEC";
 100         case GROUP_VALUE:    return "GROUP";
 101     }
 102     return NULL;
 103 }
 104
 105 /*! \copydoc _gmx_selelem_type_str() */
 106 const char *
 107 _gmx_selelem_boolean_type_str(const gmx::SelectionTreeElement &sel)
 108 {
 109     switch (sel.u.boolt)
 110     {
 111         case BOOL_NOT:  return "NOT"; break;
 112         case BOOL_AND:  return "AND"; break;
 113         case BOOL_OR:   return "OR";  break;
 114         case BOOL_XOR:  return "XOR"; break;
 115     }
 116     return NULL;
 117 }
 118
 119
 120 namespace gmx
 121 {
 122
 123 SelectionTreeElement::SelectionTreeElement(e_selelem_t type)
 124 {
 125     this->type       = type;
 126     this->flags      = (type != SEL_ROOT) ? SEL_ALLOCVAL : 0;
 127     if (type == SEL_BOOLEAN)
 128     {
 129         this->v.type = GROUP_VALUE;
 130         this->flags |= SEL_ALLOCDATA;
 131     }
 132     else
 133     {
 134         this->v.type = NO_VALUE;
 135     }
 136     _gmx_selvalue_clear(&this->v);
 137     std::memset(&this->u, 0, sizeof(this->u));
 138     this->evaluate   = NULL;
 139     this->mempool    = NULL;
 140     this->cdata      = NULL;
 141 }
 142
 143 SelectionTreeElement::~SelectionTreeElement()
 144 {
 145     /* Free the children.
 146      * Must be done before freeing other data, because the children may hold
 147      * references to data in this element. */
 148     child.reset();
 149
 150     freeValues();
 151     freeExpressionData();
 152     freeCompilerData();
 153 }
 154
 155 void SelectionTreeElement::freeValues()
 156 {
 157     mempoolRelease();
 158     if ((flags & SEL_ALLOCDATA) && v.u.ptr)
 159     {
 160         /* The number of position/group structures is constant, so the
 161          * backup of using sel->v.nr should work for them.
 162          * For strings, we report an error if we don't know the allocation
 163          * size here. */
 164         int n = (v.nalloc > 0) ? v.nalloc : v.nr;
 165         switch (v.type)
 166         {
 167             case STR_VALUE:
 168                 GMX_RELEASE_ASSERT(v.nalloc != 0,
 169                                    "SEL_ALLOCDATA should only be set for allocated "
 170                                    "STR_VALUE values");
 171                 for (int i = 0; i < n; ++i)
 172                 {
 173                     sfree(v.u.s[i]);
 174                 }
 175                 break;
 176             case GROUP_VALUE:
 177                 for (int i = 0; i < n; ++i)
 178                 {
 179                     gmx_ana_index_deinit(&v.u.g[i]);
 180                 }
 181                 break;
 182             default: /* No special handling for other types */
 183                 break;
 184         }
 185     }
 186     if (v.nalloc > 0)
 187     {
 188         if (v.type == POS_VALUE)
 189         {
 190             delete [] v.u.p;
 191         }
 192         else
 193         {
 194             sfree(v.u.ptr);
 195         }
 196     }
 197     _gmx_selvalue_setstore(&v, NULL);
 198     if (type == SEL_SUBEXPRREF && u.param)
 199     {
 200         u.param->val.u.ptr = NULL;
 201     }
 202 }
 203
 204 void
 205 SelectionTreeElement::freeExpressionData()
 206 {
 207     if (type == SEL_EXPRESSION || type == SEL_MODIFIER)
 208     {
 209         _gmx_selelem_free_method(u.expr.method, u.expr.mdata);
 210         u.expr.mdata  = NULL;
 211         u.expr.method = NULL;
 212         /* Free position data */
 213         delete u.expr.pos;
 214         u.expr.pos = NULL;
 215         /* Free position calculation data */
 216         if (u.expr.pc)
 217         {
 218             gmx_ana_poscalc_free(u.expr.pc);
 219             u.expr.pc = NULL;
 220         }
 221     }
 222     if (type == SEL_ARITHMETIC)
 223     {
 224         sfree(u.arith.opstr);
 225         u.arith.opstr = NULL;
 226     }
 227     if (type == SEL_SUBEXPR || type == SEL_ROOT
 228         || (type == SEL_CONST && v.type == GROUP_VALUE))
 229     {
 230         gmx_ana_index_deinit(&u.cgrp);
 231     }
 232     if (type == SEL_GROUPREF)
 233     {
 234         sfree(u.gref.name);
 235     }
 236 }
 237
 238 void SelectionTreeElement::mempoolReserve(int count)
 239 {
 240     if (!mempool)
 241     {
 242         return;
 243     }
 244     switch (v.type)
 245     {
 246         case INT_VALUE:
 247             v.u.i = static_cast<int *>(
 248                     _gmx_sel_mempool_alloc(mempool, sizeof(*v.u.i)*count));
 249             break;
 250
 251         case REAL_VALUE:
 252             v.u.r = static_cast<real *>(
 253                     _gmx_sel_mempool_alloc(mempool, sizeof(*v.u.r)*count));
 254             break;
 255
 256         case GROUP_VALUE:
 257             _gmx_sel_mempool_alloc_group(mempool, v.u.g, count);
 258             break;
 259
 260         default:
 261             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
 262     }
 263 }
 264
 265 void SelectionTreeElement::mempoolRelease()
 266 {
 267     if (!mempool)
 268     {
 269         return;
 270     }
 271     switch (v.type)
 272     {
 273         case INT_VALUE:
 274         case REAL_VALUE:
 275             _gmx_sel_mempool_free(mempool, v.u.ptr);
 276             _gmx_selvalue_setstore(&v, NULL);
 277             break;
 278
 279         case GROUP_VALUE:
 280             if (v.u.g)
 281             {
 282                 _gmx_sel_mempool_free_group(mempool, v.u.g);
 283             }
 284             break;
 285
 286         default:
 287             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
 288     }
 289 }
 290
 291 void SelectionTreeElement::fillNameIfMissing(const char *selectionText)
 292 {
 293     GMX_RELEASE_ASSERT(type == SEL_ROOT,
 294                        "Should not be called for non-root elements");
 295     if (name().empty())
 296     {
 297         // Check whether the actual selection given was from an external group,
 298         // and if so, use the name of the external group.
 299         SelectionTreeElementPointer child = this->child;
 300         while (child->type == SEL_MODIFIER)
 301         {
 302             if (!child->child || child->child->type != SEL_SUBEXPRREF
 303                 || !child->child->child)
 304             {
 305                 break;
 306             }
 307             child = child->child->child;
 308         }
 309         if (child->type == SEL_EXPRESSION
 310             && child->child && child->child->type == SEL_SUBEXPRREF
 311             && child->child->child)
 312         {
 313             if (child->child->child->type == SEL_CONST
 314                 && child->child->child->v.type == GROUP_VALUE)
 315             {
 316                 setName(child->child->child->name());
 317                 return;
 318             }
 319             // If the group reference is still unresolved, leave the name empty
 320             // and fill it later.
 321             if (child->child->child->type == SEL_GROUPREF)
 322             {
 323                 return;
 324             }
 325         }
 326         // If there still is no name, use the selection string.
 327         setName(selectionText);
 328     }
 329 }
 330
 331 void SelectionTreeElement::resolveIndexGroupReference(gmx_ana_indexgrps_t *grps)
 332 {
 333     GMX_RELEASE_ASSERT(type == SEL_GROUPREF,
 334                        "Should only be called for index group reference elements");
 335     if (grps == NULL)
 336     {
 337         std::string message = formatString(
 338                     "Cannot match '%s', because index groups are not available.",
 339                     name().c_str());
 340         GMX_THROW(InconsistentInputError(message));
 341     }
 342
 343     gmx_ana_index_t foundGroup;
 344     std::string     foundName;
 345     if (u.gref.name != NULL)
 346     {
 347         if (!gmx_ana_indexgrps_find(&foundGroup, &foundName, grps, u.gref.name))
 348         {
 349             std::string message = formatString(
 350                         "Cannot match '%s', because no such index group can be found.",
 351                         name().c_str());
 352             GMX_THROW(InconsistentInputError(message));
 353         }
 354     }
 355     else
 356     {
 357         if (!gmx_ana_indexgrps_extract(&foundGroup, &foundName, grps, u.gref.id))
 358         {
 359             std::string message = formatString(
 360                         "Cannot match '%s', because no such index group can be found.",
 361                         name().c_str());
 362             GMX_THROW(InconsistentInputError(message));
 363         }
 364     }
 365
 366     if (!gmx_ana_index_check_sorted(&foundGroup))
 367     {
 368         flags |= SEL_UNSORTED;
 369         // TODO: Add this test elsewhere, where it does not break valid use cases.
 370 #if 0
 371         gmx_ana_index_deinit(&foundGroup);
 372         std::string message = formatString(
 373                     "Group '%s' ('%s') cannot be used in selections, "
 374                     "because atom indices in it are not sorted and/or "
 375                     "it contains duplicate atoms.",
 376                     foundName.c_str(), name().c_str());
 377         GMX_THROW(InconsistentInputError(message));
 378 #endif
 379     }
 380
 381     sfree(u.gref.name);
 382     type = SEL_CONST;
 383     gmx_ana_index_set(&u.cgrp, foundGroup.isize, foundGroup.index,
 384                       foundGroup.nalloc_index);
 385     setName(foundName);
 386 }
 387
 388 } // namespace gmx
 389
 390 /*!
 391  * \param[in,out] sel   Selection element to set the type for.
 392  * \param[in]     vtype Value type for the selection element.
 393  *
 394  * If the new type is \ref GROUP_VALUE or \ref POS_VALUE, the
 395  * \ref SEL_ALLOCDATA flag is also set.
 396  *
 397  * This function should only be called at most once for each element,
 398  * preferably right after calling _gmx_selelem_create().
 399  */
 400 void
 401 _gmx_selelem_set_vtype(const gmx::SelectionTreeElementPointer &sel,
 402                        e_selvalue_t                            vtype)
 403 {
 404     GMX_RELEASE_ASSERT(sel->type != SEL_BOOLEAN || vtype == GROUP_VALUE,
 405                        "Boolean elements must have a group value");
 406     GMX_RELEASE_ASSERT(sel->v.type == NO_VALUE || vtype == sel->v.type,
 407                        "_gmx_selelem_set_vtype() called more than once");
 408     sel->v.type = vtype;
 409     if (vtype == GROUP_VALUE || vtype == POS_VALUE)
 410     {
 411         sel->flags |= SEL_ALLOCDATA;
 412     }
 413 }
 414
 415 /*!
 416  * \param[in] method Method to free.
 417  * \param[in] mdata  Method data to free.
 418  */
 419 void
 420 _gmx_selelem_free_method(gmx_ana_selmethod_t *method, void *mdata)
 421 {
 422     sel_freefunc free_func = NULL;
 423
 424     /* Save the pointer to the free function. */
 425     if (method && method->free)
 426     {
 427         free_func = method->free;
 428     }
 429
 430     /* Free the method itself.
 431      * Has to be done before freeing the method data, because parameter
 432      * values are typically stored in the method data, and here we may
 433      * access them. */
 434     if (method)
 435     {
 436         int  i, j;
 437
 438         /* Free the memory allocated for the parameters that are not managed
 439          * by the selection method itself. */
 440         for (i = 0; i < method->nparams; ++i)
 441         {
 442             gmx_ana_selparam_t *param = &method->param[i];
 443
 444             if (param->val.u.ptr)
 445             {
 446                 if (param->val.type == GROUP_VALUE)
 447                 {
 448                     for (j = 0; j < param->val.nr; ++j)
 449                     {
 450                         gmx_ana_index_deinit(&param->val.u.g[j]);
 451                     }
 452                 }
 453                 else if (param->val.type == POS_VALUE)
 454                 {
 455                     if (param->val.nalloc > 0)
 456                     {
 457                         delete[] param->val.u.p;
 458                         _gmx_selvalue_setstore(&param->val, NULL);
 459                     }
 460                 }
 461
 462                 if (param->val.nalloc > 0)
 463                 {
 464                     sfree(param->val.u.ptr);
 465                 }
 466             }
 467         }
 468         sfree(method->param);
 469         sfree(method);
 470     }
 471     /* Free method data. */
 472     if (mdata)
 473     {
 474         if (free_func)
 475         {
 476             free_func(mdata);
 477         }
 478         else
 479         {
 480             sfree(mdata);
 481         }
 482     }
 483 }
 484
 485 /*!
 486  * \param[in] fp      File handle to receive the output.
 487  * \param[in] sel     Root of the selection subtree to print.
 488  * \param[in] bValues If true, the evaluated values of selection elements
 489  *   are printed as well.
 490  * \param[in] level   Indentation level, starting from zero.
 491  */
 492 void
 493 _gmx_selelem_print_tree(FILE *fp, const gmx::SelectionTreeElement &sel,
 494                         bool bValues, int level)
 495 {
 496     int          i;
 497
 498     fprintf(fp, "%*c %s %s", level*2+1, '*',
 499             _gmx_selelem_type_str(sel), _gmx_sel_value_type_str(&sel.v));
 500     if (!sel.name().empty())
 501     {
 502         fprintf(fp, " \"%s\"", sel.name().c_str());
 503     }
 504     fprintf(fp, " flg=");
 505     if (sel.flags & SEL_FLAGSSET)
 506     {
 507         fprintf(fp, "s");
 508     }
 509     if (sel.flags & SEL_SINGLEVAL)
 510     {
 511         fprintf(fp, "S");
 512     }
 513     if (sel.flags & SEL_ATOMVAL)
 514     {
 515         fprintf(fp, "A");
 516     }
 517     if (sel.flags & SEL_VARNUMVAL)
 518     {
 519         fprintf(fp, "V");
 520     }
 521     if (sel.flags & SEL_DYNAMIC)
 522     {
 523         fprintf(fp, "D");
 524     }
 525     if (!(sel.flags & SEL_VALFLAGMASK))
 526     {
 527         fprintf(fp, "0");
 528     }
 529     if (sel.flags & SEL_ALLOCVAL)
 530     {
 531         fprintf(fp, "Av");
 532     }
 533     if (sel.flags & SEL_ALLOCDATA)
 534     {
 535         fprintf(fp, "Ad");
 536     }
 537     if (sel.mempool)
 538     {
 539         fprintf(fp, "P");
 540     }
 541     if (sel.type == SEL_CONST)
 542     {
 543         if (sel.v.type == INT_VALUE)
 544         {
 545             fprintf(fp, " %d", sel.v.u.i[0]);
 546         }
 547         else if (sel.v.type == REAL_VALUE)
 548         {
 549             fprintf(fp, " %f", sel.v.u.r[0]);
 550         }
 551         else if (sel.v.type == GROUP_VALUE)
 552         {
 553             const gmx_ana_index_t *g = sel.v.u.g;
 554             if (!g || g->isize == 0)
 555             {
 556                 g = &sel.u.cgrp;
 557             }
 558             fprintf(fp, " (%d atoms)", g->isize);
 559         }
 560     }
 561     else if (sel.type == SEL_BOOLEAN)
 562     {
 563         fprintf(fp, " %s", _gmx_selelem_boolean_type_str(sel));
 564     }
 565     else if (sel.type == SEL_EXPRESSION
 566              && sel.u.expr.method->name == sm_compare.name)
 567     {
 568         _gmx_selelem_print_compare_info(fp, sel.u.expr.mdata);
 569     }
 570     if (sel.evaluate)
 571     {
 572         fprintf(fp, " eval=");
 573         _gmx_sel_print_evalfunc_name(fp, sel.evaluate);
 574     }
 575     if (!(sel.flags & SEL_ALLOCVAL))
 576     {
 577         fprintf(fp, " (ext)");
 578     }
 579     fprintf(fp, "\n");
 580
 581     if ((sel.type == SEL_CONST && sel.v.type == GROUP_VALUE) || sel.type == SEL_ROOT)
 582     {
 583         const gmx_ana_index_t *g = sel.v.u.g;
 584         if (!g || g->isize == 0 || sel.evaluate != NULL)
 585         {
 586             g = &sel.u.cgrp;
 587         }
 588         if (g->isize < 0)
 589         {
 590             fprintf(fp, "%*c group: (null)\n", level*2+1, ' ');
 591         }
 592         else if (g->isize > 0)
 593         {
 594             fprintf(fp, "%*c group:", level*2+1, ' ');
 595             if (g->isize <= 20)
 596             {
 597                 for (i = 0; i < g->isize; ++i)
 598                 {
 599                     fprintf(fp, " %d", g->index[i] + 1);
 600                 }
 601             }
 602             else
 603             {
 604                 fprintf(fp, " %d atoms", g->isize);
 605             }
 606             fprintf(fp, "\n");
 607         }
 608     }
 609     else if (sel.type == SEL_EXPRESSION)
 610     {
 611         if (sel.u.expr.pc)
 612         {
 613             fprintf(fp, "%*c COM", level*2+3, '*');
 614             fprintf(fp, "\n");
 615         }
 616     }
 617
 618     if (sel.cdata)
 619     {
 620         _gmx_selelem_print_compiler_info(fp, sel, level);
 621     }
 622
 623     if (bValues && sel.type != SEL_CONST && sel.type != SEL_ROOT && sel.v.u.ptr)
 624     {
 625         fprintf(fp, "%*c value: ", level*2+1, ' ');
 626         switch (sel.v.type)
 627         {
 628             case POS_VALUE:
 629                 /* In normal use, the pointer should never be NULL, but it's
 630                  * useful to have the check for debugging to avoid accidental
 631                  * segfaults when printing the selection tree. */
 632                 if (sel.v.u.p->x)
 633                 {
 634                     fprintf(fp, "(%f, %f, %f)",
 635                             sel.v.u.p->x[0][XX], sel.v.u.p->x[0][YY],
 636                             sel.v.u.p->x[0][ZZ]);
 637                 }
 638                 else
 639                 {
 640                     fprintf(fp, "(null)");
 641                 }
 642                 break;
 643             case GROUP_VALUE:
 644                 fprintf(fp, "%d atoms", sel.v.u.g->isize);
 645                 if (sel.v.u.g->isize < 20)
 646                 {
 647                     if (sel.v.u.g->isize > 0)
 648                     {
 649                         fprintf(fp, ":");
 650                     }
 651                     for (i = 0; i < sel.v.u.g->isize; ++i)
 652                     {
 653                         fprintf(fp, " %d", sel.v.u.g->index[i] + 1);
 654                     }
 655                 }
 656                 break;
 657             default:
 658                 fprintf(fp, "???");
 659                 break;
 660         }
 661         fprintf(fp, "\n");
 662     }
 663
 664     /* Print the subexpressions with one more level of indentation */
 665     gmx::SelectionTreeElementPointer child = sel.child;
 666     while (child)
 667     {
 668         if (!(sel.type == SEL_SUBEXPRREF && child->type == SEL_SUBEXPR))
 669         {
 670             _gmx_selelem_print_tree(fp, *child, bValues, level+1);
 671         }
 672         child = child->next;
 673     }
 674 }
 675
 676 /*!
 677  * \param[in] root Root of the subtree to query.
 678  * \returns true if \p root or any any of its elements require topology
 679  *   information, false otherwise.
 680  */
 681 bool
 682 _gmx_selelem_requires_top(const gmx::SelectionTreeElement &root)
 683 {
 684     if (root.type == SEL_EXPRESSION || root.type == SEL_MODIFIER)
 685     {
 686         if (root.u.expr.method && (root.u.expr.method->flags & SMETH_REQTOP))
 687         {
 688             return true;
 689         }
 690         if (root.u.expr.pc && gmx_ana_poscalc_requires_top(root.u.expr.pc))
 691         {
 692             return true;
 693         }
 694     }
 695     gmx::SelectionTreeElementPointer child = root.child;
 696     while (child)
 697     {
 698         if (_gmx_selelem_requires_top(*child))
 699         {
 700             return true;
 701         }
 702         child = child->next;
 703     }
 704     return false;
 705 }