src/gromacs/selection/selelem.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements functions in selelem.h.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "gmx_fatal.h"
  43 #include "smalloc.h"
  44
  45 #include "gromacs/selection/indexutil.h"
  46 #include "gromacs/selection/poscalc.h"
  47 #include "gromacs/selection/position.h"
  48 #include "gromacs/selection/selmethod.h"
  49 #include "gromacs/utility/exceptions.h"
  50
  51 #include "keywords.h"
  52 #include "mempool.h"
  53 #include "selelem.h"
  54
  55 /*!
  56  * \param[in] sel Selection for which the string is requested
  57  * \returns   Pointer to a string that corresponds to \p sel->type.
  58  *
  59  * The return value points to a string constant and should not be \p free'd.
  60  *
  61  * The function returns NULL if \p sel->type is not one of the valid values.
  62  */
  63 const char *
  64 _gmx_selelem_type_str(t_selelem *sel)
  65 {
  66     switch (sel->type)
  67     {
  68         case SEL_CONST:      return "CONST";
  69         case SEL_EXPRESSION: return "EXPR";
  70         case SEL_BOOLEAN:    return "BOOL";
  71         case SEL_ARITHMETIC: return "ARITH";
  72         case SEL_ROOT:       return "ROOT";
  73         case SEL_SUBEXPR:    return "SUBEXPR";
  74         case SEL_SUBEXPRREF: return "REF";
  75         case SEL_GROUPREF:   return "GROUPREF";
  76         case SEL_MODIFIER:   return "MODIFIER";
  77     }
  78     return NULL;
  79 }
  80
  81 /*!
  82  * \param[in] val Value structore for which the string is requested.
  83  * \returns   Pointer to a string that corresponds to \p val->type,
  84  *   NULL if the type value is invalid.
  85  *
  86  * The return value points to a string constant and should not be \p free'd.
  87  */
  88 const char *
  89 _gmx_sel_value_type_str(gmx_ana_selvalue_t *val)
  90 {
  91     switch (val->type)
  92     {
  93         case NO_VALUE:       return "NONE";
  94         case INT_VALUE:      return "INT";
  95         case REAL_VALUE:     return "REAL";
  96         case STR_VALUE:      return "STR";
  97         case POS_VALUE:      return "VEC";
  98         case GROUP_VALUE:    return "GROUP";
  99     }
 100     return NULL;
 101 }
 102
 103 /*! \copydoc _gmx_selelem_type_str() */
 104 const char *
 105 _gmx_selelem_boolean_type_str(t_selelem *sel)
 106 {
 107     switch (sel->u.boolt)
 108     {
 109         case BOOL_NOT:  return "NOT"; break;
 110         case BOOL_AND:  return "AND"; break;
 111         case BOOL_OR:   return "OR";  break;
 112         case BOOL_XOR:  return "XOR"; break;
 113     }
 114     return NULL;
 115 }
 116
 117 /*!
 118  * \param[in] type Type of selection element to allocate.
 119  * \returns   Pointer to the newly allocated and initialized element.
 120  *
 121  * \c t_selelem::type is set to \p type,
 122  * \c t_selelem::v::type is set to \ref GROUP_VALUE for boolean and comparison
 123  * expressions and \ref NO_VALUE for others,
 124  * \ref SEL_ALLOCVAL is set for non-root elements (\ref SEL_ALLOCDATA is also
 125  * set for \ref SEL_BOOLEAN elements),
 126  * and \c t_selelem::refcount is set to one.
 127  * All the pointers are set to NULL.
 128  */
 129 t_selelem *
 130 _gmx_selelem_create(e_selelem_t type)
 131 {
 132     t_selelem *sel;
 133
 134     snew(sel, 1);
 135     sel->name       = NULL;
 136     sel->type       = type;
 137     sel->flags      = (type != SEL_ROOT) ? SEL_ALLOCVAL : 0;
 138     if (type == SEL_BOOLEAN)
 139     {
 140         sel->v.type = GROUP_VALUE;
 141         sel->flags |= SEL_ALLOCDATA;
 142     }
 143     else
 144     {
 145         sel->v.type = NO_VALUE;
 146     }
 147     _gmx_selvalue_clear(&sel->v);
 148     sel->evaluate   = NULL;
 149     sel->mempool    = NULL;
 150     sel->child      = NULL;
 151     sel->next       = NULL;
 152     sel->refcount   = 1;
 153
 154     return sel;
 155 }
 156
 157 /*!
 158  * \param[in,out] sel   Selection element to set the type for.
 159  * \param[in]     vtype Value type for the selection element.
 160  * \returns       0 on success, EINVAL if the value type is invalid.
 161  *
 162  * If the new type is \ref GROUP_VALUE or \ref POS_VALUE, the
 163  * \ref SEL_ALLOCDATA flag is also set.
 164  *
 165  * This function should only be called at most once for each element,
 166  * preferably right after calling _gmx_selelem_create().
 167  */
 168 int
 169 _gmx_selelem_set_vtype(t_selelem *sel, e_selvalue_t vtype)
 170 {
 171     if (sel->type == SEL_BOOLEAN && vtype != GROUP_VALUE)
 172     {
 173         gmx_bug("internal error");
 174         return EINVAL;
 175     }
 176     if (sel->v.type != NO_VALUE && vtype != sel->v.type)
 177     {
 178         gmx_call("_gmx_selelem_set_vtype() called more than once");
 179         return EINVAL;
 180     }
 181     sel->v.type = vtype;
 182     if (vtype == GROUP_VALUE || vtype == POS_VALUE)
 183     {
 184         sel->flags |= SEL_ALLOCDATA;
 185     }
 186     return 0;
 187 }
 188
 189 /*!
 190  * \param[in,out] sel   Selection element to reserve.
 191  * \param[in]     count Number of values to reserve memory for.
 192  * \returns       0 on success or if no memory pool, non-zero on error.
 193  *
 194  * Reserves memory for the values of \p sel from the \p sel->mempool
 195  * memory pool. If no memory pool is set, nothing is done.
 196  */
 197 void
 198 _gmx_selelem_mempool_reserve(t_selelem *sel, int count)
 199 {
 200     if (!sel->mempool)
 201     {
 202         return;
 203     }
 204     switch (sel->v.type)
 205     {
 206         case INT_VALUE:
 207             sel->v.u.i = static_cast<int *>(
 208                     _gmx_sel_mempool_alloc(sel->mempool, sizeof(*sel->v.u.i)*count));
 209             break;
 210
 211         case REAL_VALUE:
 212             sel->v.u.r = static_cast<real *>(
 213                     _gmx_sel_mempool_alloc(sel->mempool, sizeof(*sel->v.u.r)*count));
 214             break;
 215
 216         case GROUP_VALUE:
 217             _gmx_sel_mempool_alloc_group(sel->mempool, sel->v.u.g, count);
 218             break;
 219
 220         default:
 221             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
 222     }
 223 }
 224
 225 /*!
 226  * \param[in,out] sel   Selection element to release.
 227  *
 228  * Releases the memory allocated for the values of \p sel from the
 229  * \p sel->mempool memory pool. If no memory pool is set, nothing is done.
 230  */
 231 void
 232 _gmx_selelem_mempool_release(t_selelem *sel)
 233 {
 234     if (!sel->mempool)
 235     {
 236         return;
 237     }
 238     switch (sel->v.type)
 239     {
 240         case INT_VALUE:
 241         case REAL_VALUE:
 242             _gmx_sel_mempool_free(sel->mempool, sel->v.u.ptr);
 243             _gmx_selvalue_setstore(&sel->v, NULL);
 244             break;
 245
 246         case GROUP_VALUE:
 247             if (sel->v.u.g)
 248             {
 249                 _gmx_sel_mempool_free_group(sel->mempool, sel->v.u.g);
 250             }
 251             break;
 252
 253         default:
 254             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
 255     }
 256 }
 257
 258 /*!
 259  * \param[in] sel Selection to free.
 260  */
 261 void
 262 _gmx_selelem_free_values(t_selelem *sel)
 263 {
 264     int   i, n;
 265
 266     _gmx_selelem_mempool_release(sel);
 267     if ((sel->flags & SEL_ALLOCDATA) && sel->v.u.ptr)
 268     {
 269         /* The number of position/group structures is constant, so the
 270          * backup of using sel->v.nr should work for them.
 271          * For strings, we report an error if we don't know the allocation
 272          * size here. */
 273         n = (sel->v.nalloc > 0) ? sel->v.nalloc : sel->v.nr;
 274         switch (sel->v.type)
 275         {
 276             case STR_VALUE:
 277                 if (sel->v.nalloc == 0)
 278                 {
 279                     gmx_bug("SEL_ALLOCDATA should only be set for allocated STR_VALUE values");
 280                     break;
 281                 }
 282                 for (i = 0; i < n; ++i)
 283                 {
 284                     sfree(sel->v.u.s[i]);
 285                 }
 286                 break;
 287             case POS_VALUE:
 288                 for (i = 0; i < n; ++i)
 289                 {
 290                     gmx_ana_pos_deinit(&sel->v.u.p[i]);
 291                 }
 292                 break;
 293             case GROUP_VALUE:
 294                 for (i = 0; i < n; ++i)
 295                 {
 296                     gmx_ana_index_deinit(&sel->v.u.g[i]);
 297                 }
 298                 break;
 299             default: /* No special handling for other types */
 300                 break;
 301         }
 302     }
 303     if (sel->flags & SEL_ALLOCVAL)
 304     {
 305         sfree(sel->v.u.ptr);
 306     }
 307     _gmx_selvalue_setstore(&sel->v, NULL);
 308     if (sel->type == SEL_SUBEXPRREF && sel->u.param)
 309     {
 310         sel->u.param->val.u.ptr = NULL;
 311     }
 312 }
 313
 314 /*!
 315  * \param[in] method Method to free.
 316  * \param[in] mdata  Method data to free.
 317  */
 318 void
 319 _gmx_selelem_free_method(gmx_ana_selmethod_t *method, void *mdata)
 320 {
 321     sel_freefunc free_func = NULL;
 322
 323     /* Save the pointer to the free function. */
 324     if (method && method->free)
 325     {
 326         free_func = method->free;
 327     }
 328
 329     /* Free the method itself.
 330      * Has to be done before freeing the method data, because parameter
 331      * values are typically stored in the method data, and here we may
 332      * access them. */
 333     if (method)
 334     {
 335         int  i, j;
 336
 337         /* Free the memory allocated for the parameters that are not managed
 338          * by the selection method itself. */
 339         for (i = 0; i < method->nparams; ++i)
 340         {
 341             gmx_ana_selparam_t *param = &method->param[i];
 342
 343             if (param->val.u.ptr)
 344             {
 345                 if (param->val.type == GROUP_VALUE)
 346                 {
 347                     for (j = 0; j < param->val.nr; ++j)
 348                     {
 349                         gmx_ana_index_deinit(&param->val.u.g[j]);
 350                     }
 351                 }
 352                 else if (param->val.type == POS_VALUE)
 353                 {
 354                     for (j = 0; j < param->val.nr; ++j)
 355                     {
 356                         gmx_ana_pos_deinit(&param->val.u.p[j]);
 357                     }
 358                 }
 359
 360                 if (param->val.nalloc > 0)
 361                 {
 362                     sfree(param->val.u.ptr);
 363                 }
 364             }
 365         }
 366         sfree(method->param);
 367         sfree(method);
 368     }
 369     /* Free method data. */
 370     if (mdata)
 371     {
 372         if (free_func)
 373         {
 374             free_func(mdata);
 375         }
 376         sfree(mdata);
 377     }
 378 }
 379
 380 /*!
 381  * \param[in] sel Selection to free.
 382  */
 383 void
 384 _gmx_selelem_free_exprdata(t_selelem *sel)
 385 {
 386     if (sel->type == SEL_EXPRESSION || sel->type == SEL_MODIFIER)
 387     {
 388         _gmx_selelem_free_method(sel->u.expr.method, sel->u.expr.mdata);
 389         sel->u.expr.mdata = NULL;
 390         sel->u.expr.method = NULL;
 391         /* Free position data */
 392         if (sel->u.expr.pos)
 393         {
 394             gmx_ana_pos_free(sel->u.expr.pos);
 395             sel->u.expr.pos = NULL;
 396         }
 397         /* Free position calculation data */
 398         if (sel->u.expr.pc)
 399         {
 400             gmx_ana_poscalc_free(sel->u.expr.pc);
 401             sel->u.expr.pc = NULL;
 402         }
 403     }
 404     if (sel->type == SEL_ARITHMETIC)
 405     {
 406         sfree(sel->u.arith.opstr);
 407         sel->u.arith.opstr = NULL;
 408     }
 409     if (sel->type == SEL_SUBEXPR || sel->type == SEL_ROOT
 410         || (sel->type == SEL_CONST && sel->v.type == GROUP_VALUE))
 411     {
 412         gmx_ana_index_deinit(&sel->u.cgrp);
 413     }
 414     if (sel->type == SEL_GROUPREF)
 415     {
 416         sfree(sel->u.gref.name);
 417     }
 418 }
 419
 420 /*!
 421  * \param[in] sel Selection to free.
 422  *
 423  * Decrements \ref t_selelem::refcount "sel->refcount" and frees the
 424  * memory allocated for \p sel and all its children if the reference count
 425  * reaches zero.
 426  */
 427 void
 428 _gmx_selelem_free(t_selelem *sel)
 429 {
 430     /* Decrement the reference counter and do nothing if references remain */
 431     sel->refcount--;
 432     if (sel->refcount > 0)
 433     {
 434         return;
 435     }
 436
 437     /* Free the children.
 438      * Must be done before freeing other data, because the children may hold
 439      * references to data in this element. */
 440     _gmx_selelem_free_chain(sel->child);
 441
 442     /* Free value storage */
 443     _gmx_selelem_free_values(sel);
 444
 445     /* Free other storage */
 446     _gmx_selelem_free_exprdata(sel);
 447
 448     /* Free temporary compiler data if present */
 449     _gmx_selelem_free_compiler_data(sel);
 450
 451     sfree(sel);
 452 }
 453
 454 /*!
 455  * \param[in] first First selection to free.
 456  *
 457  * Frees \p first and all selections accessible through the
 458  * \ref t_selelem::next "first->next" pointer.
 459  */
 460 void
 461 _gmx_selelem_free_chain(t_selelem *first)
 462 {
 463     t_selelem *child, *prev;
 464
 465     child = first;
 466     while (child)
 467     {
 468         prev = child;
 469         child = child->next;
 470         _gmx_selelem_free(prev);
 471     }
 472 }
 473
 474 /*!
 475  * \param[in] fp      File handle to receive the output.
 476  * \param[in] sel     Root of the selection subtree to print.
 477  * \param[in] bValues If true, the evaluated values of selection elements
 478  *   are printed as well.
 479  * \param[in] level   Indentation level, starting from zero.
 480  */
 481 void
 482 _gmx_selelem_print_tree(FILE *fp, t_selelem *sel, bool bValues, int level)
 483 {
 484     t_selelem *child;
 485     int          i;
 486
 487     fprintf(fp, "%*c %s %s", level*2+1, '*',
 488             _gmx_selelem_type_str(sel), _gmx_sel_value_type_str(&sel->v));
 489     if (sel->name)
 490     {
 491         fprintf(fp, " \"%s\"", sel->name);
 492     }
 493     fprintf(fp, " flg=");
 494     if (sel->flags & SEL_FLAGSSET)
 495     {
 496         fprintf(fp, "s");
 497     }
 498     if (sel->flags & SEL_SINGLEVAL)
 499     {
 500         fprintf(fp, "S");
 501     }
 502     if (sel->flags & SEL_ATOMVAL)
 503     {
 504         fprintf(fp, "A");
 505     }
 506     if (sel->flags & SEL_VARNUMVAL)
 507     {
 508         fprintf(fp, "V");
 509     }
 510     if (sel->flags & SEL_DYNAMIC)
 511     {
 512         fprintf(fp, "D");
 513     }
 514     if (!(sel->flags & SEL_VALFLAGMASK))
 515     {
 516         fprintf(fp, "0");
 517     }
 518     if (sel->mempool)
 519     {
 520         fprintf(fp, "P");
 521     }
 522     if (sel->type == SEL_CONST)
 523     {
 524         if (sel->v.type == INT_VALUE)
 525         {
 526             fprintf(fp, " %d", sel->v.u.i[0]);
 527         }
 528         else if (sel->v.type == REAL_VALUE)
 529         {
 530             fprintf(fp, " %f", sel->v.u.r[0]);
 531         }
 532         else if (sel->v.type == GROUP_VALUE)
 533         {
 534             gmx_ana_index_t *g = sel->v.u.g;
 535             if (!g || g->isize == 0)
 536                 g = &sel->u.cgrp;
 537             fprintf(fp, " (%d atoms)", g->isize);
 538         }
 539     }
 540     else if (sel->type == SEL_BOOLEAN)
 541     {
 542         fprintf(fp, " %s", _gmx_selelem_boolean_type_str(sel));
 543     }
 544     else if (sel->type == SEL_EXPRESSION
 545              && sel->u.expr.method->name == sm_compare.name)
 546     {
 547         _gmx_selelem_print_compare_info(fp, sel->u.expr.mdata);
 548     }
 549     if (sel->evaluate)
 550     {
 551         fprintf(fp, " eval=");
 552         _gmx_sel_print_evalfunc_name(fp, sel->evaluate);
 553     }
 554     if (sel->refcount > 1)
 555     {
 556         fprintf(fp, " refc=%d", sel->refcount);
 557     }
 558     if (!(sel->flags & SEL_ALLOCVAL))
 559     {
 560         fprintf(fp, " (ext. output)");
 561     }
 562     fprintf(fp, "\n");
 563
 564     if ((sel->type == SEL_CONST && sel->v.type == GROUP_VALUE) || sel->type == SEL_ROOT)
 565     {
 566         gmx_ana_index_t *g = sel->v.u.g;
 567         if (!g || g->isize == 0 || sel->evaluate != NULL)
 568         {
 569             g = &sel->u.cgrp;
 570         }
 571         if (g->isize < 0)
 572         {
 573             fprintf(fp, "%*c group: (null)\n", level*2+1, ' ');
 574         }
 575         else if (g->isize > 0)
 576         {
 577             fprintf(fp, "%*c group:", level*2+1, ' ');
 578             if (g->isize <= 20)
 579             {
 580                 for (i = 0; i < g->isize; ++i)
 581                 {
 582                     fprintf(fp, " %d", g->index[i] + 1);
 583                 }
 584             }
 585             else
 586             {
 587                 fprintf(fp, " %d atoms", g->isize);
 588             }
 589             fprintf(fp, "\n");
 590         }
 591     }
 592     else if (sel->type == SEL_EXPRESSION)
 593     {
 594         if (sel->u.expr.pc)
 595         {
 596             fprintf(fp, "%*c COM", level*2+3, '*');
 597             fprintf(fp, "\n");
 598         }
 599     }
 600
 601     if (sel->cdata)
 602     {
 603         _gmx_selelem_print_compiler_info(fp, sel, level);
 604     }
 605
 606     if (bValues && sel->type != SEL_CONST && sel->type != SEL_ROOT && sel->v.u.ptr)
 607     {
 608         fprintf(fp, "%*c value: ", level*2+1, ' ');
 609         switch (sel->v.type)
 610         {
 611             case POS_VALUE:
 612                 /* In normal use, the pointer should never be NULL, but it's
 613                  * useful to have the check for debugging to avoid accidental
 614                  * segfaults when printing the selection tree. */
 615                 if (sel->v.u.p->x)
 616                 {
 617                     fprintf(fp, "(%f, %f, %f)",
 618                             sel->v.u.p->x[0][XX], sel->v.u.p->x[0][YY],
 619                             sel->v.u.p->x[0][ZZ]);
 620                 }
 621                 else
 622                 {
 623                     fprintf(fp, "(null)");
 624                 }
 625                 break;
 626             case GROUP_VALUE:
 627                 fprintf(fp, "%d atoms", sel->v.u.g->isize);
 628                 if (sel->v.u.g->isize < 20)
 629                 {
 630                     if (sel->v.u.g->isize > 0)
 631                     {
 632                         fprintf(fp, ":");
 633                     }
 634                     for (i = 0; i < sel->v.u.g->isize; ++i)
 635                     {
 636                         fprintf(fp, " %d", sel->v.u.g->index[i] + 1);
 637                     }
 638                 }
 639                 break;
 640             default:
 641                 fprintf(fp, "???");
 642                 break;
 643         }
 644         fprintf(fp, "\n");
 645     }
 646
 647     /* Print the subexpressions with one more level of indentation */
 648     child = sel->child;
 649     while (child)
 650     {
 651         if (!(sel->type == SEL_SUBEXPRREF && child->type == SEL_SUBEXPR))
 652         {
 653             _gmx_selelem_print_tree(fp, child, bValues, level+1);
 654         }
 655         child = child->next;
 656     }
 657 }
 658
 659 /*!
 660  * \param[in] root Root of the subtree to query.
 661  * \returns true if \p root or any any of its elements require topology
 662  *   information, false otherwise.
 663  */
 664 bool
 665 _gmx_selelem_requires_top(t_selelem *root)
 666 {
 667     t_selelem *child;
 668
 669     if (root->type == SEL_EXPRESSION || root->type == SEL_MODIFIER)
 670     {
 671         if (root->u.expr.method && (root->u.expr.method->flags & SMETH_REQTOP))
 672         {
 673             return true;
 674         }
 675         if (root->u.expr.pc && gmx_ana_poscalc_requires_top(root->u.expr.pc))
 676         {
 677             return true;
 678         }
 679     }
 680     child = root->child;
 681     while (child)
 682     {
 683         if (_gmx_selelem_requires_top(child))
 684         {
 685             return true;
 686         }
 687         child = child->next;
 688     }
 689     return false;
 690 }