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