src/gromacs/selection/compiler.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 Selection compilation and optimization.
  33  *
  34  * \todo
  35  * Better error handling and memory management in error situations.
  36  * At least, the main compilation function leaves the selection collection in
  37  * a bad state if an error occurs.
  38  *
  39  * \todo
  40  * The memory usage could still be optimized.
  41  * Use of memory pooling could still be extended, and a lot of redundant
  42  * gmin/gmax data could be eliminated for complex arithmetic expressions.
  43  *
  44  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  45  * \ingroup module_selection
  46  */
  47 /*! \internal
  48  * \page page_module_selection_compiler Selection compilation
  49  *
  50  * The compiler takes the selection element tree from the selection parser
  51  * (see \ref page_module_selection_parser) as input.
  52  * The selection parser is quite independent of selection evaluation details,
  53  * and the compiler processes the tree to conform to what the evaluation
  54  * functions expect.
  55  * For better control and optimization possibilities, the compilation is
  56  * done on all selections simultaneously.
  57  * Hence, all the selections should be parsed before the compiler can be
  58  * called.
  59  *
  60  * The compiler initializes all fields in gmx::SelectionTreeElement not
  61  * initialized by the parser: gmx::SelectionTreeElement::v (some fields have
  62  * already been initialized by the parser),
  63  * gmx::SelectionTreeElement::evaluate, and gmx::SelectionTreeElement::u
  64  * (again, some elements have been initialized in the parser).
  65  * The gmx::SelectionTreeElement::cdata field is used during the compilation to store
  66  * internal data, but the data is freed when the compiler returns.
  67  *
  68  * In addition to initializing the elements, the compiler reorganizes the tree
  69  * to simplify and optimize evaluation. The compiler also evaluates the static
  70  * parts of the selection: in the end of the compilation, static parts have
  71  * been replaced by the result of the evaluation.
  72  *
  73  * The compiler is invoked using gmx::SelectionCompiler.
  74  * The gmx::SelectionCompiler::compile() method does the compilation in several
  75  * passes over the gmx::SelectionTreeElement tree.
  76  *  -# Defaults are set for the position type and flags of position calculation
  77  *     methods that were not explicitly specified in the user input.
  78  *  -# Subexpressions are extracted: a separate root is created for each
  79  *     subexpression, and placed before the expression is first used.
  80  *     Currently, only variables and expressions used to evaluate parameter
  81  *     values are extracted, but common subexpression could also be detected
  82  *     here.
  83  *  -# A second pass (in fact, multiple passes because of interdependencies)
  84  *     with simple reordering and initialization is done:
  85  *    -# Boolean expressions are combined such that one element can evaluate,
  86  *       e.g., "A and B and C". The subexpressions in boolean expression are
  87  *       reordered such that static expressions come first without otherwise
  88  *       altering the relative order of the expressions.
  89  *    -# The compiler data structure is allocated for each element, and
  90  *       the fields are initialized, with the exception of the contents of
  91  *       \c gmax and \c gmin fields.  This is the part that needs multiple
  92  *       passes, because some flags are set recursively based on which elements
  93  *       refer to an element, and these flags need to be set to initialize
  94  *       other fields.
  95  *    -# The gmx::SelectionTreeElement::evaluate field is set to the correct
  96  *       evaluation function from evaluate.h.
  97  *    .
  98  *  -# The evaluation function of all elements is replaced with the
  99  *     analyze_static() function to be able to initialize the element before
 100  *     the actual evaluation function is called.
 101  *     The evaluation machinery is then called to initialize the whole tree,
 102  *     while simultaneously evaluating the static expressions.
 103  *     During the evaluation, track is kept of the smallest and largest
 104  *     possible selections, and these are stored in the internal compiler
 105  *     data structure for each element.
 106  *     To be able to do this for all possible values of dynamical expressions,
 107  *     special care needs to be taken with boolean expressions because they
 108  *     are short-circuiting. This is done through the
 109  *     \c SEL_CDATA_EVALMAX flag, which makes dynamic child expressions
 110  *     of \c BOOL_OR expressions evaluate to empty groups, while subexpressions
 111  *     of \c BOOL_AND are evaluated to largest possible groups.
 112  *     Memory is also allocated to store the results of the evaluation.
 113  *     For each element, analyze_static() calls the actual evaluation function
 114  *     after the element has been properly initialized.
 115  *  -# Another evaluation pass is done over subexpressions with more than
 116  *     one reference to them. These cannot be completely processed during the
 117  *     first pass, because it is not known whether later references require
 118  *     additional evaluation of static expressions.
 119  *  -# Unused subexpressions are removed. For efficiency reasons (and to avoid
 120  *     some checks), this is actually done several times already earlier in
 121  *     the compilation process.
 122  *  -# Most of the processing is now done, and the next pass simply sets the
 123  *     evaluation group of root elements to the largest selection as determined
 124  *     in pass 4.  For root elements of subexpressions that should not be
 125  *     evaluated before they are referred to, the evaluation group/function is
 126  *     cleared.  At the same time, position calculation data is initialized for
 127  *     for selection method elements that require it.  Compiler data is also
 128  *     freed as it is no longer needed.
 129  *  -# A final pass initializes the total masses and charges in the
 130  *     \c gmx_ana_selection_t data structures.
 131  *
 132  * The actual evaluation of the selection is described in the documentation
 133  * of the functions in evaluate.h.
 134  *
 135  * \todo
 136  * Some combinations of method parameter flags are not yet properly treated by
 137  * the compiler or the evaluation functions in evaluate.cpp. All the ones used by
 138  * currently implemented methods should work, but new combinations might not.
 139  *
 140  *
 141  * \section selcompiler_tree Element tree after compilation
 142  *
 143  * After the compilation, the selection element tree is suitable for
 144  * gmx_ana_selcollection_evaluate().
 145  * Enough memory has been allocated for gmx::SelectionTreeElement::v
 146  * (and gmx::SelectionTreeElement::cgrp for \ref SEL_SUBEXPR elements) to allow
 147  * the selection to be evaluated without allocating any memory.
 148  *
 149  *
 150  * \subsection selcompiler_tree_root Root elements
 151  *
 152  * The top level of the tree consists of a chain of \ref SEL_ROOT elements.
 153  * These are used for two purposes:
 154  *  -# A selection that should be evaluated.
 155  *     These elements appear in the same order as the selections in the input.
 156  *     For these elements, gmx::SelectionTreeElement::v has been set to the
 157  *     maximum possible group that the selection can evaluate to (only for
 158  *     dynamic selections), and gmx::SelectionTreeElement::cgrp has been set to
 159  *     use a NULL group for evaluation.
 160  *  -# A subexpression that appears in one or more selections.
 161  *     Each selection that gives a value for a method parameter is a
 162  *     potential subexpression, as is any variable value.
 163  *     Only subexpressions that require evaluation for each frame are left
 164  *     after the selection is compiled.
 165  *     Each subexpression appears in the chain before any references to it.
 166  *     For these elements, gmx::SelectionTreeElement::cgrp has been set to the
 167  *     group that should be used to evaluate the subexpression.
 168  *     If gmx::SelectionTreeElement::cgrp is empty, the total evaluation group
 169  *     is not known in advance or it is more efficient to evaluate the
 170  *     subexpression only when it is referenced.  If this is the case,
 171  *     gmx::SelectionTreeElement::evaluate is also NULL.
 172  *
 173  * The children of the \ref SEL_ROOT elements can be used to distinguish
 174  * the two types of root elements from each other; the rules are the same
 175  * as for the parsed tree (see \ref selparser_tree_root).
 176  * Subexpressions are treated as if they had been provided through variables.
 177  *
 178  * Selection names are stored as after parsing (see \ref selparser_tree_root).
 179  *
 180  *
 181  * \subsection selcompiler_tree_const Constant elements
 182  *
 183  * All (sub)selections that do not require particle positions have been
 184  * replaced with \ref SEL_CONST elements.
 185  * Constant elements from the parser are also retained if present in
 186  * dynamic parts of the selections.
 187  * Several constant elements with a NULL gmx::SelectionTreeElement::evaluate
 188  * are left for debugging purposes; of these, only the ones for \ref BOOL_OR
 189  * expressions are used during evaluation.
 190  *
 191  * The value is stored in gmx::SelectionTreeElement::v, and for group values
 192  * with an evaluation function set, also in gmx::SelectionTreeElement::cgrp.
 193  * For \ref GROUP_VALUE elements, unnecessary atoms (i.e., atoms that
 194  * could never be selected) have been removed from the value.
 195  *
 196  * \ref SEL_CONST elements have no children.
 197  *
 198  *
 199  * \subsection selcompiler_tree_method Method evaluation elements
 200  *
 201  * All selection methods that need to be evaluated dynamically are described
 202  * by a \ref SEL_EXPRESSION element. The gmx::SelectionTreeElement::method and
 203  * gmx::SelectionTreeElement::mdata fields have already been initialized by the parser,
 204  * and the compiler only calls the initialization functions in the method
 205  * data structure to do some additional initialization of these fields at
 206  * appropriate points. If the gmx::SelectionTreeElement::pc data field has been
 207  * created by the parser, the compiler initializes the data structure properly
 208  * once the required positions are known. If the gmx::SelectionTreeElement::pc
 209  * field is NULL after the parser, but the method provides only
 210  * sel_updatefunc_pos(), an appropriate position calculation data structure is
 211  * created.  If gmx::SelectionTreeElement::pc is not NULL,
 212  * gmx::SelectionTreeElement::pos is also initialized to hold the positions
 213  * calculated.
 214  *
 215  * Children of these elements are of type \ref SEL_SUBEXPRREF, and describe
 216  * parameter values that need to be evaluated for each frame. See the next
 217  * section for more details.
 218  * \ref SEL_CONST children can also appear, and stand for parameters that get
 219  * their value from a static expression. These elements are present only for
 220  * debugging purposes: they always have a NULL evaluation function.
 221  *
 222  *
 223  * \subsection selcompiler_tree_subexpr Subexpression elements
 224  *
 225  * As described in \ref selcompiler_tree_root, subexpressions are created
 226  * for each variable and each expression that gives a value to a selection
 227  * method parameter. As the only child of the \ref SEL_ROOT element,
 228  * these elements have a \ref SEL_SUBEXPR element. The \ref SEL_SUBEXPR
 229  * element has a single child, which evaluates the actual expression.
 230  * After compilation, only subexpressions that require particle positions
 231  * for evaluation are left.
 232  * For non-variable subexpression, automatic names have been generated to
 233  * help in debugging.
 234  *
 235  * For \ref SEL_SUBEXPR elements, memory has been allocated for
 236  * gmx::SelectionTreeElement::cgrp to store the group for which the expression
 237  * has been evaluated during the current frame.  This is only done if full
 238  * subexpression evaluation by _gmx_sel_evaluate_subexpr() is needed; the other
 239  * evaluation functions do not require this memory.
 240  *
 241  * \ref SEL_SUBEXPRREF elements are used to describe references to
 242  * subexpressions. They have always a single child, which is the
 243  * \ref SEL_SUBEXPR element being referenced.
 244  *
 245  * If a subexpression is used only once, the evaluation has been optimized by
 246  * setting the child of the \ref SEL_SUBEXPR element to evaluate the value of
 247  * \ref SEL_SUBEXPRREF directly (in the case of memory pooling, this is managed
 248  * by the evaluation functions).  In such cases, the evaluation routines for the
 249  * \ref SEL_SUBEXPRREF and \ref SEL_SUBEXPR elements only propagate some status
 250  * information, but do not unnecessarily copy the values.
 251  *
 252  *
 253  * \subsection selcompiler_tree_bool Boolean elements
 254  *
 255  * \ref SEL_BOOLEAN elements have been merged such that one element
 256  * may carry out evaluation of more than one operation of the same type.
 257  * The static parts of the expressions have been evaluated, and are placed
 258  * in the first child. These are followed by the dynamic expressions, in the
 259  * order provided by the user.
 260  *
 261  *
 262  * \subsection selcompiler_tree_arith Arithmetic elements
 263  *
 264  * Constant and static expressions in \ref SEL_ARITHMETIC elements have been
 265  * calculated.
 266  * Currently, no other processing is done.
 267  */
 268 #include "compiler.h"
 269
 270 #ifdef HAVE_CONFIG_H
 271 #include <config.h>
 272 #endif
 273
 274 #include <algorithm>
 275
 276 #include <math.h>
 277 #include <stdarg.h>
 278
 279 #include "smalloc.h"
 280 #include "string2.h"
 281 #include "vec.h"
 282
 283 #include "gromacs/selection/indexutil.h"
 284 #include "gromacs/selection/poscalc.h"
 285 #include "gromacs/selection/selection.h"
 286 #include "gromacs/selection/selmethod.h"
 287 #include "gromacs/utility/exceptions.h"
 288 #include "gromacs/utility/stringutil.h"
 289
 290 #include "evaluate.h"
 291 #include "keywords.h"
 292 #include "mempool.h"
 293 #include "selectioncollection-impl.h"
 294 #include "selelem.h"
 295
 296 using std::min;
 297 using gmx::SelectionTreeElement;
 298 using gmx::SelectionTreeElementPointer;
 299
 300 /*! \internal \brief
 301  * Compiler flags.
 302  */
 303 enum
 304 {
 305     /*! \brief
 306      * Whether a subexpression needs to evaluated for all atoms.
 307      *
 308      * This flag is set for \ref SEL_SUBEXPR elements that are used to
 309      * evaluate non-atom-valued selection method parameters, as well as
 310      * those that are used directly as values of selections.
 311      */
 312     SEL_CDATA_FULLEVAL    =  1,
 313     /*! \brief
 314      * Whether the whole subexpression should be treated as static.
 315      *
 316      * This flag is always false if \ref SEL_DYNAMIC is set for the element,
 317      * but it is also false for static elements within common subexpressions.
 318      */
 319     SEL_CDATA_STATIC      =  2,
 320     /** Whether the subexpression will always be evaluated in the same group. */
 321     SEL_CDATA_STATICEVAL  =  4,
 322     /** Whether the compiler evaluation routine should return the maximal selection. */
 323     SEL_CDATA_EVALMAX     =  8,
 324     /** Whether memory has been allocated for \p gmin and \p gmax. */
 325     SEL_CDATA_MINMAXALLOC = 16,
 326     /** Whether to update \p gmin and \p gmax in static analysis. */
 327     SEL_CDATA_DOMINMAX      = 128,
 328     /** Whether the subexpression uses simple pass evaluation functions. */
 329     SEL_CDATA_SIMPLESUBEXPR = 32,
 330     /** Whether this expression is a part of a common subexpression. */
 331     SEL_CDATA_COMMONSUBEXPR = 64
 332 };
 333
 334 /*! \internal \brief
 335  * Internal data structure used by the compiler.
 336  */
 337 typedef struct t_compiler_data
 338 {
 339     /** The real evaluation method. */
 340     gmx::sel_evalfunc evaluate;
 341     /** Number of references to a \ref SEL_SUBEXPR element. */
 342     int              refcount;
 343     /** Flags for specifying how to treat this element during compilation. */
 344     int              flags;
 345     /** Smallest selection that can be selected by the subexpression. */
 346     gmx_ana_index_t *gmin;
 347     /** Largest selection that can be selected by the subexpression. */
 348     gmx_ana_index_t *gmax;
 349 } t_compiler_data;
 350
 351
 352 /********************************************************************
 353  * COMPILER UTILITY FUNCTIONS
 354  ********************************************************************/
 355
 356 /*! \brief
 357  * Helper method for printing out debug information about a min/max group.
 358  */
 359 static void
 360 print_group_info(FILE *fp, const char *name,
 361                  const SelectionTreeElement &sel,
 362                  gmx_ana_index_t *g)
 363 {
 364     fprintf(fp, " %s=", name);
 365     if (!g)
 366     {
 367         fprintf(fp, "(null)");
 368     }
 369     else if (sel.cdata->flags & SEL_CDATA_MINMAXALLOC)
 370     {
 371         fprintf(fp, "(%d atoms, %p)", g->isize, (void*)g);
 372     }
 373     else if (sel.v.type == GROUP_VALUE && g == sel.v.u.g)
 374     {
 375         fprintf(fp, "(static, %p)", (void*)g);
 376     }
 377     else
 378     {
 379         fprintf(fp, "%p", (void*)g);
 380     }
 381 }
 382
 383 /*!
 384  * \param[in] fp      File handle to receive the output.
 385  * \param[in] sel     Selection element to print.
 386  * \param[in] level   Indentation level, starting from zero.
 387  */
 388 void
 389 _gmx_selelem_print_compiler_info(FILE *fp, const SelectionTreeElement &sel,
 390                                  int level)
 391 {
 392     if (!sel.cdata)
 393     {
 394         return;
 395     }
 396     fprintf(fp, "%*c cdata: flg=", level*2+1, ' ');
 397     if (sel.cdata->flags & SEL_CDATA_FULLEVAL)
 398     {
 399         fprintf(fp, "F");
 400     }
 401     if (!(sel.cdata->flags & SEL_CDATA_STATIC))
 402     {
 403         fprintf(fp, "D");
 404     }
 405     if (sel.cdata->flags & SEL_CDATA_STATICEVAL)
 406     {
 407         fprintf(fp, "S");
 408     }
 409     if (sel.cdata->flags & SEL_CDATA_EVALMAX)
 410     {
 411         fprintf(fp, "M");
 412     }
 413     if (sel.cdata->flags & SEL_CDATA_MINMAXALLOC)
 414     {
 415         fprintf(fp, "A");
 416     }
 417     if (sel.cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
 418     {
 419         fprintf(fp, "Ss");
 420     }
 421     if (sel.cdata->flags & SEL_CDATA_COMMONSUBEXPR)
 422     {
 423         fprintf(fp, "Sc");
 424     }
 425     if (!sel.cdata->flags)
 426     {
 427         fprintf(fp, "0");
 428     }
 429     if (sel.cdata->refcount > 0)
 430     {
 431         fprintf(fp, " refc=%d", sel.cdata->refcount);
 432     }
 433     fprintf(fp, " eval=");
 434     _gmx_sel_print_evalfunc_name(fp, sel.cdata->evaluate);
 435     print_group_info(fp, "gmin", sel, sel.cdata->gmin);
 436     print_group_info(fp, "gmax", sel, sel.cdata->gmax);
 437     fprintf(fp, "\n");
 438 }
 439
 440 namespace gmx
 441 {
 442
 443 void SelectionTreeElement::freeCompilerData()
 444 {
 445     if (cdata)
 446     {
 447         evaluate = cdata->evaluate;
 448         if (cdata->flags & SEL_CDATA_MINMAXALLOC)
 449         {
 450             cdata->gmin->name = NULL;
 451             cdata->gmax->name = NULL;
 452             gmx_ana_index_deinit(cdata->gmin);
 453             gmx_ana_index_deinit(cdata->gmax);
 454             sfree(cdata->gmin);
 455             sfree(cdata->gmax);
 456         }
 457         sfree(cdata);
 458     }
 459     cdata = NULL;
 460 }
 461
 462 } // namespace gmx
 463
 464 /*! \brief
 465  * Allocates memory for storing the evaluated value of a selection element.
 466  *
 467  * \param     sel   Selection element to initialize
 468  * \param[in] isize Maximum evaluation group size.
 469  * \param[in] bChildEval true if children have already been processed.
 470  *
 471  * If called more than once, memory is (re)allocated to ensure that the
 472  * maximum of the \p isize values can be stored.
 473  */
 474 static void
 475 alloc_selection_data(const SelectionTreeElementPointer &sel,
 476                      int isize, bool bChildEval)
 477 {
 478     int        nalloc;
 479
 480     if (sel->mempool)
 481     {
 482         return;
 483     }
 484     /* Find out the number of elements to allocate */
 485     if (sel->flags & SEL_SINGLEVAL)
 486     {
 487         nalloc = 1;
 488     }
 489     else if (sel->flags & SEL_ATOMVAL)
 490     {
 491         nalloc = isize;
 492     }
 493     else /* sel->flags should contain SEL_VARNUMVAL */
 494     {
 495         if (!bChildEval)
 496         {
 497             return;
 498         }
 499         SelectionTreeElementPointer child = sel;
 500         if (sel->type == SEL_SUBEXPRREF)
 501         {
 502             GMX_RELEASE_ASSERT(sel->child && sel->child->type == SEL_SUBEXPR,
 503                 "Subexpression expected for subexpression reference");
 504             child = sel->child->child;
 505             GMX_RELEASE_ASSERT(child,
 506                 "Subexpression elements should always have a child element");
 507         }
 508         nalloc = (sel->v.type == POS_VALUE) ? child->v.u.p->nr : child->v.nr;
 509     }
 510     /* For positions, we actually want to allocate just a single structure
 511      * for nalloc positions. */
 512     if (sel->v.type == POS_VALUE)
 513     {
 514         isize  = nalloc;
 515         nalloc = 1;
 516     }
 517     /* Allocate memory for sel->v.u if needed */
 518     if (sel->flags & SEL_ALLOCVAL)
 519     {
 520         _gmx_selvalue_reserve(&sel->v, nalloc);
 521     }
 522     /* Reserve memory inside group and position structures if
 523      * SEL_ALLOCDATA is set. */
 524     if (sel->flags & SEL_ALLOCDATA)
 525     {
 526         if (sel->v.type == GROUP_VALUE)
 527         {
 528             gmx_ana_index_reserve(sel->v.u.g, isize);
 529         }
 530         else if (sel->v.type == POS_VALUE)
 531         {
 532             gmx_ana_pos_reserve(sel->v.u.p, isize, 0);
 533         }
 534     }
 535 }
 536
 537 /*! \brief
 538  * Replace the evaluation function of each element in the subtree.
 539  *
 540  * \param     sel  Root of the selection subtree to process.
 541  * \param[in] eval The new evaluation function.
 542  */
 543 static void
 544 set_evaluation_function(const SelectionTreeElementPointer &sel,
 545                         gmx::sel_evalfunc eval)
 546 {
 547     sel->evaluate = eval;
 548     if (sel->type != SEL_SUBEXPRREF)
 549     {
 550         SelectionTreeElementPointer child = sel->child;
 551         while (child)
 552         {
 553             set_evaluation_function(child, eval);
 554             child = child->next;
 555         }
 556     }
 557 }
 558
 559
 560 /********************************************************************
 561  * POSITION KEYWORD DEFAULT INITIALIZATION
 562  ********************************************************************/
 563
 564 /*! \brief
 565  * Initializes default values for position keyword evaluation.
 566  *
 567  * \param[in,out] root       Root of the element tree to initialize.
 568  * \param[in]     spost      Default output position type.
 569  * \param[in]     rpost      Default reference position type.
 570  * \param[in]     sel        Selection that the element evaluates the positions
 571  *      for, or NULL if the element is an internal element.
 572  */
 573 static void
 574 init_pos_keyword_defaults(SelectionTreeElement *root,
 575                           const char *spost, const char *rpost,
 576                           const gmx::internal::SelectionData *sel)
 577 {
 578     /* Selections use largest static group by default, while
 579      * reference positions use the whole residue/molecule. */
 580     if (root->type == SEL_EXPRESSION)
 581     {
 582         bool bSelection = (sel != NULL);
 583         int flags = bSelection ? POS_COMPLMAX : POS_COMPLWHOLE;
 584         if (bSelection)
 585         {
 586             if (sel->hasFlag(gmx::efSelection_DynamicMask))
 587             {
 588                 flags |= POS_MASKONLY;
 589             }
 590             if (sel->hasFlag(gmx::efSelection_EvaluateVelocities))
 591             {
 592                 flags |= POS_VELOCITIES;
 593             }
 594             if (sel->hasFlag(gmx::efSelection_EvaluateForces))
 595             {
 596                 flags |= POS_FORCES;
 597             }
 598         }
 599         _gmx_selelem_set_kwpos_type(root, bSelection ? spost : rpost);
 600         _gmx_selelem_set_kwpos_flags(root, flags);
 601     }
 602     /* Change the defaults once we are no longer processing modifiers */
 603     if (root->type != SEL_ROOT && root->type != SEL_MODIFIER
 604         && root->type != SEL_SUBEXPRREF && root->type != SEL_SUBEXPR)
 605     {
 606         sel = NULL;
 607     }
 608     /* Recurse into children */
 609     SelectionTreeElementPointer child = root->child;
 610     while (child)
 611     {
 612         init_pos_keyword_defaults(child.get(), spost, rpost, sel);
 613         child = child->next;
 614     }
 615 }
 616
 617
 618 /********************************************************************
 619  * SUBEXPRESSION PROCESSING
 620  ********************************************************************/
 621
 622 /*! \brief
 623  * Reverses the chain of selection elements starting at \p root.
 624  *
 625  * \param   root First selection in the whole selection chain.
 626  * \returns The new first element for the chain.
 627  */
 628 static SelectionTreeElementPointer
 629 reverse_selelem_chain(const SelectionTreeElementPointer &root)
 630 {
 631     SelectionTreeElementPointer prev;
 632     SelectionTreeElementPointer item = root;
 633     while (item)
 634     {
 635         SelectionTreeElementPointer next = item->next;
 636         item->next = prev;
 637         prev = item;
 638         item = next;
 639     }
 640     return prev;
 641 }
 642
 643 /*! \brief
 644  * Removes subexpressions that don't have any references.
 645  *
 646  * \param     root First selection in the whole selection chain.
 647  * \returns   The new first element for the chain.
 648  *
 649  * The elements are processed in reverse order to correctly detect
 650  * subexpressions only referred to by other subexpressions.
 651  */
 652 static SelectionTreeElementPointer
 653 remove_unused_subexpressions(SelectionTreeElementPointer root)
 654 {
 655     if (!root)
 656     {
 657         return SelectionTreeElementPointer();
 658     }
 659     root = reverse_selelem_chain(root);
 660     while (root->child->type == SEL_SUBEXPR && root->child.unique())
 661     {
 662         // Frees the root element.
 663         root = root->next;
 664     }
 665     SelectionTreeElementPointer prev = root;
 666     SelectionTreeElementPointer item = root->next;
 667     while (item)
 668     {
 669         SelectionTreeElementPointer next = item->next;
 670         if (item->child->type == SEL_SUBEXPR && item->child.unique())
 671         {
 672             // Frees the current item when it goes out of scope.
 673             prev->next = next;
 674         }
 675         else
 676         {
 677             prev = item;
 678         }
 679         item = next;
 680     }
 681     return reverse_selelem_chain(root);
 682 }
 683
 684 /*! \brief
 685  * Creates a name with a running number for a subexpression.
 686  *
 687  * \param[in,out] sel The subexpression to be named.
 688  * \param[in]     i   Running number for the subexpression.
 689  *
 690  * The name of the selection becomes "SubExpr N", where N is \p i;
 691  * Memory is allocated for the name and the name is stored both in
 692  * gmx::SelectionTreeElement::name and gmx::SelectionTreeElement::u::cgrp::name; the latter
 693  * is freed by _gmx_selelem_free().
 694  */
 695 static void
 696 create_subexpression_name(const SelectionTreeElementPointer &sel, int i)
 697 {
 698     char *name = strdup(gmx::formatString("SubExpr %d", i).c_str());
 699
 700     sel->name        = name;
 701     sel->u.cgrp.name = name;
 702 }
 703
 704 /*! \brief
 705  * Processes and extracts subexpressions from a given selection subtree.
 706  *
 707  * \param   sel      Root of the subtree to process.
 708  * \param   subexprn Pointer to a subexpression counter.
 709  * \returns Pointer to a chain of subselections, or NULL if none were found.
 710  *
 711  * This function finds recursively all \ref SEL_SUBEXPRREF elements below
 712  * the given root element and ensures that their children are within
 713  * \ref SEL_SUBEXPR elements. It also creates a chain of \ref SEL_ROOT elements
 714  * that contain the subexpression as their children and returns the first
 715  * of these root elements.
 716  */
 717 static SelectionTreeElementPointer
 718 extract_item_subselections(const SelectionTreeElementPointer &sel,
 719                            int *subexprn)
 720 {
 721     SelectionTreeElementPointer root;
 722     SelectionTreeElementPointer subexpr;
 723     SelectionTreeElementPointer child = sel->child;
 724
 725     while (child)
 726     {
 727         if (!root)
 728         {
 729             root = subexpr = extract_item_subselections(child, subexprn);
 730         }
 731         else
 732         {
 733             subexpr->next = extract_item_subselections(child, subexprn);
 734         }
 735         while (subexpr && subexpr->next)
 736         {
 737             subexpr = subexpr->next;
 738         }
 739         /* The latter check excludes variable references.
 740          * It also excludes subexpression elements that have already been
 741          * processed, because they are given a name when they are first
 742          * encountered.
 743          * TODO: There should be a more robust mechanism (probably a dedicated
 744          * flag) for detecting parser-generated subexpressions than relying on
 745          * a NULL name field. Additional TODO: This mechanism doesn't seem to
 746          * be currently used. */
 747         if (child->type == SEL_SUBEXPRREF && (child->child->type != SEL_SUBEXPR
 748                                               || child->child->name == NULL))
 749         {
 750             /* Create the root element for the subexpression */
 751             if (!root)
 752             {
 753                 root.reset(new SelectionTreeElement(SEL_ROOT));
 754                 subexpr = root;
 755             }
 756             else
 757             {
 758                 subexpr->next.reset(new SelectionTreeElement(SEL_ROOT));
 759                 subexpr = subexpr->next;
 760             }
 761             /* Create the subexpression element and/or
 762              * move the actual subexpression under the created element. */
 763             if (child->child->type != SEL_SUBEXPR)
 764             {
 765                 subexpr->child.reset(new SelectionTreeElement(SEL_SUBEXPR));
 766                 _gmx_selelem_set_vtype(subexpr->child, child->v.type);
 767                 subexpr->child->child = child->child;
 768                 child->child          = subexpr->child;
 769             }
 770             else
 771             {
 772                 subexpr->child = child->child;
 773             }
 774             create_subexpression_name(subexpr->child, ++*subexprn);
 775             /* Set the flags for the created elements */
 776             subexpr->flags          |= (child->flags & SEL_VALFLAGMASK);
 777             subexpr->child->flags   |= (child->flags & SEL_VALFLAGMASK);
 778         }
 779         child = child->next;
 780     }
 781
 782     return root;
 783 }
 784
 785 /*! \brief
 786  * Extracts subexpressions of the selection chain.
 787  *
 788  * \param   sel First selection in the whole selection chain.
 789  * \returns The new first element for the chain.
 790  *
 791  * Finds all the subexpressions (and their subexpressions) in the
 792  * selection chain starting from \p sel and creates \ref SEL_SUBEXPR
 793  * elements for them.
 794  * \ref SEL_ROOT elements are also created for each subexpression
 795  * and inserted into the selection chain before the expressions that
 796  * refer to them.
 797  */
 798 static SelectionTreeElementPointer
 799 extract_subexpressions(SelectionTreeElementPointer sel)
 800 {
 801     SelectionTreeElementPointer root;
 802     SelectionTreeElementPointer next = sel;
 803     int subexprn = 0;
 804     while (next)
 805     {
 806         SelectionTreeElementPointer item
 807             = extract_item_subselections(next, &subexprn);
 808         if (item)
 809         {
 810             if (!root)
 811             {
 812                 root = item;
 813             }
 814             else
 815             {
 816                 sel->next = item;
 817             }
 818             while (item->next)
 819             {
 820                 item = item->next;
 821             }
 822             item->next = next;
 823         }
 824         else if (!root)
 825         {
 826             root = next;
 827         }
 828         sel = next;
 829         next = next->next;
 830     }
 831     return root;
 832 }
 833
 834
 835 /********************************************************************
 836  * BOOLEAN OPERATION REORDERING
 837  ********************************************************************/
 838
 839 /*! \brief
 840  * Removes redundant boolean selection elements.
 841  *
 842  * \param  sel Root of the selection subtree to optimize.
 843  *
 844  * This function merges similar boolean operations (e.g., (A or B) or C becomes
 845  * a single OR operation with three operands).
 846  */
 847 static void
 848 optimize_boolean_expressions(const SelectionTreeElementPointer &sel)
 849 {
 850     /* Do recursively for children */
 851     if (sel->type != SEL_SUBEXPRREF)
 852     {
 853         SelectionTreeElementPointer prev;
 854         SelectionTreeElementPointer child = sel->child;
 855         while (child)
 856         {
 857             optimize_boolean_expressions(child);
 858             /* Remove double negations */
 859             if (child->type == SEL_BOOLEAN && child->u.boolt == BOOL_NOT
 860                 && child->child->type == SEL_BOOLEAN && child->child->u.boolt == BOOL_NOT)
 861             {
 862                 /* Move the doubly negated expression up two levels */
 863                 if (!prev)
 864                 {
 865                     sel->child = child->child->child;
 866                     prev       = sel->child;
 867                 }
 868                 else
 869                 {
 870                     prev->next = child->child->child;
 871                     prev       = prev->next;
 872                 }
 873                 child->child->child->next = child->next;
 874                 // Discards the two negations.
 875                 child = prev;
 876             }
 877             prev  = child;
 878             child = child->next;
 879         }
 880     }
 881     if (sel->type != SEL_BOOLEAN || sel->u.boolt == BOOL_NOT)
 882     {
 883         return;
 884     }
 885     /* Merge subsequent binary operations */
 886     SelectionTreeElementPointer prev;
 887     SelectionTreeElementPointer child = sel->child;
 888     while (child)
 889     {
 890         if (child->type == SEL_BOOLEAN && child->u.boolt == sel->u.boolt)
 891         {
 892             if (!prev)
 893             {
 894                 sel->child = child->child;
 895                 prev       = sel->child;
 896             }
 897             else
 898             {
 899                 prev->next = child->child;
 900             }
 901             while (prev->next)
 902             {
 903                 prev = prev->next;
 904             }
 905             prev->next = child->next;
 906             // Discards the old child.
 907             child = prev->next;
 908         }
 909         else
 910         {
 911             prev = child;
 912             child = child->next;
 913         }
 914     }
 915 }
 916
 917 /*! \brief
 918  * Reorders children of boolean expressions such that static selections
 919  * come first.
 920  *
 921  * \param  sel Root of the selection subtree to reorder.
 922  *
 923  * The relative order of static expressions does not change.
 924  * The same is true for the dynamic expressions.
 925  */
 926 static void
 927 reorder_boolean_static_children(const SelectionTreeElementPointer &sel)
 928 {
 929     /* Do recursively for children */
 930     if (sel->type != SEL_SUBEXPRREF)
 931     {
 932         SelectionTreeElementPointer child = sel->child;
 933         while (child)
 934         {
 935             reorder_boolean_static_children(child);
 936             child = child->next;
 937         }
 938     }
 939
 940     /* Reorder boolean expressions such that static selections come first */
 941     if (sel->type == SEL_BOOLEAN && (sel->flags & SEL_DYNAMIC))
 942     {
 943         // Add a dummy head element that precedes the first child.
 944         SelectionTreeElementPointer dummy(
 945                 new SelectionTreeElement(SEL_BOOLEAN));
 946         dummy->next = sel->child;
 947         SelectionTreeElementPointer prev  = dummy;
 948         SelectionTreeElementPointer child = dummy;
 949         while (child->next)
 950         {
 951             /* child is the last handled static expression */
 952             /* prev is the last handled non-static expression */
 953             SelectionTreeElementPointer next = prev->next;
 954             while (next && (next->flags & SEL_DYNAMIC))
 955             {
 956                 prev = next;
 957                 next = next->next;
 958             }
 959             /* next is now the first static expression after child */
 960             if (!next)
 961             {
 962                 break;
 963             }
 964             /* Reorder such that next comes after child */
 965             if (prev != child)
 966             {
 967                 prev->next  = next->next;
 968                 next->next  = child->next;
 969                 child->next = next;
 970             }
 971             else
 972             {
 973                 prev = prev->next;
 974             }
 975             /* Advance child by one */
 976             child = next;
 977         }
 978
 979         sel->child = dummy->next;
 980     }
 981 }
 982
 983
 984 /********************************************************************
 985  * ARITHMETIC EXPRESSION PROCESSING
 986  ********************************************************************/
 987
 988 /*! \brief
 989  * Processes arithmetic expressions to simplify and speed up evaluation.
 990  *
 991  * \param  sel Root of the selection subtree to process.
 992  *
 993  * Currently, this function only converts integer constants to reals
 994  * within arithmetic expressions.
 995  */
 996 static void
 997 optimize_arithmetic_expressions(const SelectionTreeElementPointer &sel)
 998 {
 999     /* Do recursively for children. */
1000     if (sel->type != SEL_SUBEXPRREF)
1001     {
1002         SelectionTreeElementPointer child = sel->child;
1003         while (child)
1004         {
1005             optimize_arithmetic_expressions(child);
1006             child = child->next;
1007         }
1008     }
1009
1010     if (sel->type != SEL_ARITHMETIC)
1011     {
1012         return;
1013     }
1014
1015     /* Convert integer constants to reals. */
1016     SelectionTreeElementPointer child = sel->child;
1017     while (child)
1018     {
1019         if (child->v.type == INT_VALUE)
1020         {
1021             real  *r;
1022
1023             if (child->type != SEL_CONST)
1024             {
1025                 GMX_THROW(gmx::InconsistentInputError("Non-constant integer expressions not implemented in arithmetic evaluation"));
1026             }
1027             snew(r, 1);
1028             r[0] = child->v.u.i[0];
1029             sfree(child->v.u.i);
1030             child->v.u.r = r;
1031             child->v.type = REAL_VALUE;
1032         }
1033         else if (child->v.type != REAL_VALUE)
1034         {
1035             GMX_THROW(gmx::InternalError("Non-numerical value in arithmetic expression"));
1036         }
1037         child = child->next;
1038     }
1039 }
1040
1041
1042 /********************************************************************
1043  * EVALUATION PREPARATION COMPILER
1044  ********************************************************************/
1045
1046 /*! \brief
1047  * Sets the evaluation functions for the selection (sub)tree.
1048  *
1049  * \param[in,out] sel Root of the selection subtree to process.
1050  *
1051  * This function sets the evaluation function
1052  * (gmx::SelectionTreeElement::evaluate) for the selection elements.
1053  */
1054 static void
1055 init_item_evalfunc(const SelectionTreeElementPointer &sel)
1056 {
1057     /* Process children. */
1058     if (sel->type != SEL_SUBEXPRREF)
1059     {
1060         SelectionTreeElementPointer child = sel->child;
1061         while (child)
1062         {
1063             init_item_evalfunc(child);
1064             child = child->next;
1065         }
1066     }
1067
1068     /* Set the evaluation function */
1069     switch (sel->type)
1070     {
1071         case SEL_CONST:
1072             if (sel->v.type == GROUP_VALUE)
1073             {
1074                 sel->evaluate = &_gmx_sel_evaluate_static;
1075             }
1076             break;
1077
1078         case SEL_EXPRESSION:
1079             if (!(sel->flags & SEL_DYNAMIC) && sel->u.expr.method
1080                 && sel->u.expr.method->init_frame)
1081             {
1082                 sel->flags |= SEL_INITFRAME;
1083             }
1084             sel->evaluate = &_gmx_sel_evaluate_method;
1085             break;
1086
1087         case SEL_ARITHMETIC:
1088             sel->evaluate = &_gmx_sel_evaluate_arithmetic;
1089             break;
1090
1091         case SEL_MODIFIER:
1092             if (sel->v.type != NO_VALUE)
1093             {
1094                 sel->evaluate = &_gmx_sel_evaluate_modifier;
1095             }
1096             break;
1097
1098         case SEL_BOOLEAN:
1099             switch (sel->u.boolt)
1100             {
1101                 case BOOL_NOT: sel->evaluate = &_gmx_sel_evaluate_not; break;
1102                 case BOOL_AND: sel->evaluate = &_gmx_sel_evaluate_and; break;
1103                 case BOOL_OR:  sel->evaluate = &_gmx_sel_evaluate_or;  break;
1104                 case BOOL_XOR:
1105                     GMX_THROW(gmx::NotImplementedError("xor expressions not implemented"));
1106             }
1107             break;
1108
1109         case SEL_ROOT:
1110             sel->evaluate = &_gmx_sel_evaluate_root;
1111             break;
1112
1113         case SEL_SUBEXPR:
1114             sel->evaluate = ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
1115                              ? &_gmx_sel_evaluate_subexpr_simple
1116                              : &_gmx_sel_evaluate_subexpr);
1117             break;
1118
1119         case SEL_SUBEXPRREF:
1120             sel->name     = sel->child->name;
1121             sel->evaluate = ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
1122                              ? &_gmx_sel_evaluate_subexprref_simple
1123                              : &_gmx_sel_evaluate_subexprref);
1124             break;
1125
1126         case SEL_GROUPREF:
1127             GMX_THROW(gmx::APIError("Unresolved group reference in compilation"));
1128     }
1129     sel->cdata->evaluate = sel->evaluate;
1130 }
1131
1132 /*! \brief
1133  * Sets the memory pool for selection elements that can use it.
1134  *
1135  * \param     sel      Root of the selection subtree to process.
1136  * \param[in] mempool  Memory pool to use.
1137  */
1138 static void
1139 setup_memory_pooling(const SelectionTreeElementPointer &sel,
1140                      gmx_sel_mempool_t *mempool)
1141 {
1142     if (sel->type != SEL_SUBEXPRREF)
1143     {
1144         SelectionTreeElementPointer child = sel->child;
1145         while (child)
1146         {
1147             if ((sel->type == SEL_BOOLEAN && (child->flags & SEL_DYNAMIC))
1148                 || (sel->type == SEL_ARITHMETIC && child->type != SEL_CONST
1149                     && !(child->flags & SEL_SINGLEVAL))
1150                 || (sel->type == SEL_SUBEXPR
1151                     && !(sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)))
1152             {
1153                 child->mempool = mempool;
1154                 if (child->type == SEL_SUBEXPRREF
1155                     && (child->child->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
1156                 {
1157                     child->child->child->mempool = mempool;
1158                 }
1159             }
1160             setup_memory_pooling(child, mempool);
1161             child = child->next;
1162         }
1163     }
1164 }
1165
1166 /*! \brief
1167  * Prepares the selection (sub)tree for evaluation.
1168  *
1169  * \param[in,out] sel Root of the selection subtree to prepare.
1170  *
1171  * It also allocates memory for the \p sel->v.u.g or \p sel->v.u.p
1172  * structure if required.
1173  */
1174 static void
1175 init_item_evaloutput(const SelectionTreeElementPointer &sel)
1176 {
1177     GMX_ASSERT(!(sel->child == NULL &&
1178                  (sel->type == SEL_SUBEXPRREF || sel->type == SEL_SUBEXPR)),
1179                "Subexpression elements should always have a child element");
1180
1181     /* Process children. */
1182     if (sel->type != SEL_SUBEXPRREF)
1183     {
1184         SelectionTreeElementPointer child = sel->child;
1185         while (child)
1186         {
1187             init_item_evaloutput(child);
1188             child = child->next;
1189         }
1190     }
1191
1192     if (sel->type == SEL_SUBEXPR
1193         && (sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
1194     {
1195         sel->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
1196         if (sel->v.type == GROUP_VALUE || sel->v.type == POS_VALUE)
1197         {
1198             _gmx_selvalue_setstore(&sel->v, sel->child->v.u.ptr);
1199         }
1200     }
1201     else if (sel->type == SEL_SUBEXPR
1202              && (sel->cdata->flags & SEL_CDATA_FULLEVAL))
1203     {
1204         sel->evaluate = &_gmx_sel_evaluate_subexpr_staticeval;
1205         sel->cdata->evaluate = sel->evaluate;
1206         sel->child->mempool = NULL;
1207         sel->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
1208         if (sel->v.type == GROUP_VALUE || sel->v.type == POS_VALUE)
1209         {
1210             _gmx_selvalue_setstore(&sel->v, sel->child->v.u.ptr);
1211         }
1212     }
1213     else if (sel->type == SEL_SUBEXPRREF
1214              && (sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
1215     {
1216         if (sel->v.u.ptr)
1217         {
1218             _gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
1219             sel->child->child->freeValues();
1220             sel->child->child->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
1221             sel->child->child->flags |= (sel->flags & SEL_ALLOCDATA);
1222             _gmx_selvalue_setstore(&sel->child->child->v, sel->v.u.ptr);
1223         }
1224         else if (sel->v.type == GROUP_VALUE || sel->v.type == POS_VALUE)
1225         {
1226             _gmx_selvalue_setstore(&sel->v, sel->child->child->v.u.ptr);
1227         }
1228         sel->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
1229     }
1230
1231     /* Make sure that the group/position structure is allocated. */
1232     if (!sel->v.u.ptr && (sel->flags & SEL_ALLOCVAL))
1233     {
1234         if (sel->v.type == GROUP_VALUE || sel->v.type == POS_VALUE)
1235         {
1236             _gmx_selvalue_reserve(&sel->v, 1);
1237             sel->v.nr = 1;
1238         }
1239     }
1240 }
1241
1242
1243 /********************************************************************
1244  * COMPILER DATA INITIALIZATION
1245  ********************************************************************/
1246
1247 /*! \brief
1248  * Allocates memory for the compiler data and initializes the structure.
1249  *
1250  * \param sel Root of the selection subtree to process.
1251  */
1252 static void
1253 init_item_compilerdata(const SelectionTreeElementPointer &sel)
1254 {
1255     /* Allocate the compiler data structure */
1256     snew(sel->cdata, 1);
1257
1258     /* Store the real evaluation method because the compiler will replace it */
1259     sel->cdata->evaluate = sel->evaluate;
1260
1261     /* This will be computed separately. */
1262     sel->cdata->refcount = 0;
1263
1264     /* Initialize the flags */
1265     sel->cdata->flags = SEL_CDATA_STATICEVAL;
1266     if (!(sel->flags & SEL_DYNAMIC))
1267     {
1268         sel->cdata->flags |= SEL_CDATA_STATIC;
1269     }
1270     if (sel->type == SEL_SUBEXPR)
1271     {
1272         sel->cdata->flags |= SEL_CDATA_EVALMAX;
1273     }
1274     /* Set the full evaluation flag for subexpressions that require it;
1275      * the subexpression has already been initialized, so we can simply
1276      * access its compilation flags.*/
1277     if (sel->type == SEL_EXPRESSION || sel->type == SEL_MODIFIER)
1278     {
1279         SelectionTreeElementPointer child = sel->child;
1280         while (child)
1281         {
1282             if (!(child->flags & SEL_ATOMVAL) && child->child)
1283             {
1284                 child->child->cdata->flags |= SEL_CDATA_FULLEVAL;
1285             }
1286             child = child->next;
1287         }
1288     }
1289     else if (sel->type == SEL_ROOT && sel->child->type == SEL_SUBEXPRREF)
1290     {
1291         sel->child->child->cdata->flags |= SEL_CDATA_FULLEVAL;
1292     }
1293
1294     /* Initialize children */
1295     if (sel->type != SEL_SUBEXPRREF)
1296     {
1297         SelectionTreeElementPointer child = sel->child;
1298         while (child)
1299         {
1300             init_item_compilerdata(child);
1301             child = child->next;
1302         }
1303     }
1304
1305     /* Determine whether we should evaluate the minimum or the maximum
1306      * for the children of this element. */
1307     if (sel->type == SEL_BOOLEAN)
1308     {
1309         bool  bEvalMax;
1310
1311         bEvalMax = (sel->u.boolt == BOOL_AND);
1312         SelectionTreeElementPointer child = sel->child;
1313         while (child)
1314         {
1315             if (bEvalMax)
1316             {
1317                 child->cdata->flags |= SEL_CDATA_EVALMAX;
1318             }
1319             else if (child->type == SEL_BOOLEAN && child->u.boolt == BOOL_NOT)
1320             {
1321                 child->child->cdata->flags |= SEL_CDATA_EVALMAX;
1322             }
1323             child = child->next;
1324         }
1325     }
1326     else if (sel->type == SEL_EXPRESSION || sel->type == SEL_MODIFIER
1327              || sel->type == SEL_SUBEXPR)
1328     {
1329         SelectionTreeElementPointer child = sel->child;
1330         while (child)
1331         {
1332             child->cdata->flags |= SEL_CDATA_EVALMAX;
1333             child = child->next;
1334         }
1335     }
1336 }
1337
1338 /*! \brief
1339  * Initializes the static evaluation flag for a selection subtree.
1340  *
1341  * \param[in,out] sel  Root of the selection subtree to process.
1342  *
1343  * Sets the \c bStaticEval in the compiler data structure:
1344  * for any element for which the evaluation group may depend on the trajectory
1345  * frame, the flag is cleared.
1346  *
1347  * reorder_boolean_static_children() should have been called.
1348  */
1349 static void
1350 init_item_staticeval(const SelectionTreeElementPointer &sel)
1351 {
1352     /* Subexpressions with full evaluation should always have bStaticEval,
1353      * so don't do anything if a reference to them is encountered. */
1354     if (sel->type == SEL_SUBEXPRREF
1355         && (sel->child->cdata->flags & SEL_CDATA_FULLEVAL))
1356     {
1357         return;
1358     }
1359
1360     /* Propagate the bStaticEval flag to children if it is not set */
1361     if (!(sel->cdata->flags & SEL_CDATA_STATICEVAL))
1362     {
1363         SelectionTreeElementPointer child = sel->child;
1364         while (child)
1365         {
1366             if ((sel->type != SEL_EXPRESSION && sel->type != SEL_MODIFIER)
1367                 || (child->flags & SEL_ATOMVAL))
1368             {
1369                 if (child->cdata->flags & SEL_CDATA_STATICEVAL)
1370                 {
1371                     child->cdata->flags &= ~SEL_CDATA_STATICEVAL;
1372                     init_item_staticeval(child);
1373                 }
1374             }
1375             child = child->next;
1376         }
1377     }
1378     else /* bStaticEval is set */
1379     {
1380         /* For boolean expressions, any expression after the first dynamic
1381          * expression should not have bStaticEval. */
1382         if (sel->type == SEL_BOOLEAN)
1383         {
1384             SelectionTreeElementPointer child = sel->child;
1385             while (child && !(child->flags & SEL_DYNAMIC))
1386             {
1387                 child = child->next;
1388             }
1389             if (child)
1390             {
1391                 child = child->next;
1392             }
1393             while (child)
1394             {
1395                 child->cdata->flags &= ~SEL_CDATA_STATICEVAL;
1396                 child = child->next;
1397             }
1398         }
1399
1400         /* Process the children */
1401         SelectionTreeElementPointer child = sel->child;
1402         while (child)
1403         {
1404             init_item_staticeval(child);
1405             child = child->next;
1406         }
1407     }
1408 }
1409
1410 /*! \brief
1411  * Compute reference counts for subexpressions.
1412  *
1413  * \param sel Root of the selection subtree to process.
1414  */
1415 static void
1416 init_item_subexpr_refcount(const SelectionTreeElementPointer &sel)
1417 {
1418     // Reset the counter when the subexpression is first encountered.
1419     if (sel->type == SEL_ROOT && sel->child->type == SEL_SUBEXPR
1420         && sel->child->cdata)
1421     {
1422         sel->child->cdata->refcount = 0;
1423     }
1424
1425     if (sel->type == SEL_SUBEXPRREF)
1426     {
1427         ++sel->child->cdata->refcount;
1428     }
1429     else
1430     {
1431         SelectionTreeElementPointer child = sel->child;
1432         while (child)
1433         {
1434             init_item_subexpr_refcount(child);
1435             child = child->next;
1436         }
1437     }
1438 }
1439
1440 /*! \brief
1441  * Initializes compiler flags for subexpressions.
1442  *
1443  * \param sel Root of the selection subtree to process.
1444  */
1445 static void
1446 init_item_subexpr_flags(const SelectionTreeElementPointer &sel)
1447 {
1448     if (sel->type == SEL_SUBEXPR)
1449     {
1450         if (sel->cdata->refcount == 1)
1451         {
1452             sel->cdata->flags |= SEL_CDATA_SIMPLESUBEXPR;
1453         }
1454         else if (!(sel->cdata->flags & SEL_CDATA_FULLEVAL))
1455         {
1456             sel->cdata->flags |= SEL_CDATA_COMMONSUBEXPR;
1457         }
1458     }
1459     else if (sel->type == SEL_SUBEXPRREF
1460              && (sel->child->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
1461     {
1462         sel->cdata->flags |= SEL_CDATA_SIMPLESUBEXPR;
1463     }
1464
1465     /* Process children, but only follow subexpression references if the
1466      * common subexpression flag needs to be propagated. */
1467     if (sel->type != SEL_SUBEXPRREF
1468         || ((sel->cdata->flags & SEL_CDATA_COMMONSUBEXPR)
1469             && !(sel->child->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)))
1470     {
1471         SelectionTreeElementPointer child = sel->child;
1472         while (child)
1473         {
1474             if (!(child->cdata->flags & SEL_CDATA_COMMONSUBEXPR))
1475             {
1476                 if (sel->type != SEL_EXPRESSION || (child->flags & SEL_ATOMVAL))
1477                 {
1478                     child->cdata->flags |=
1479                         (sel->cdata->flags & SEL_CDATA_COMMONSUBEXPR);
1480                 }
1481                 init_item_subexpr_flags(child);
1482             }
1483             child = child->next;
1484         }
1485     }
1486 }
1487
1488 /*! \brief
1489  * Initializes the gmin and gmax fields of the compiler data structure.
1490  *
1491  * \param sel Root of the selection subtree to process.
1492  */
1493 static void
1494 init_item_minmax_groups(const SelectionTreeElementPointer &sel)
1495 {
1496     /* Process children. */
1497     if (sel->type != SEL_SUBEXPRREF)
1498     {
1499         SelectionTreeElementPointer child = sel->child;
1500         while (child)
1501         {
1502             init_item_minmax_groups(child);
1503             child = child->next;
1504         }
1505     }
1506
1507     /* Initialize the minimum and maximum evaluation groups. */
1508     if (sel->type != SEL_ROOT && sel->v.type != NO_VALUE)
1509     {
1510         if (sel->v.type == GROUP_VALUE
1511             && (sel->cdata->flags & SEL_CDATA_STATIC))
1512         {
1513             sel->cdata->gmin = sel->v.u.g;
1514             sel->cdata->gmax = sel->v.u.g;
1515         }
1516         else if (sel->type == SEL_SUBEXPR
1517                  && ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
1518                      || (sel->cdata->flags & SEL_CDATA_FULLEVAL)))
1519         {
1520             GMX_ASSERT(sel->child,
1521                        "Subexpression elements should always have a child element");
1522             sel->cdata->gmin = sel->child->cdata->gmin;
1523             sel->cdata->gmax = sel->child->cdata->gmax;
1524         }
1525         else
1526         {
1527             sel->cdata->flags |= SEL_CDATA_MINMAXALLOC | SEL_CDATA_DOMINMAX;
1528             snew(sel->cdata->gmin, 1);
1529             snew(sel->cdata->gmax, 1);
1530         }
1531     }
1532 }
1533
1534
1535 /********************************************************************
1536  * EVALUATION GROUP INITIALIZATION
1537  ********************************************************************/
1538
1539 /*! \brief
1540  * Initializes evaluation groups for root items.
1541  *
1542  * \param[in,out] sc   Selection collection data.
1543  *
1544  * The evaluation group of each \ref SEL_ROOT element corresponding to a
1545  * selection in \p sc is set to \p gall.  The same is done for \ref SEL_ROOT
1546  * elements corresponding to subexpressions that need full evaluation.
1547  */
1548 static void
1549 initialize_evalgrps(gmx_ana_selcollection_t *sc)
1550 {
1551     SelectionTreeElementPointer root = sc->root;
1552     while (root)
1553     {
1554         GMX_RELEASE_ASSERT(root->child,
1555                            "Root elements should always have a child");
1556         if (root->child->type != SEL_SUBEXPR
1557             || (root->child->cdata->flags & SEL_CDATA_FULLEVAL))
1558         {
1559             gmx_ana_index_set(&root->u.cgrp, sc->gall.isize, sc->gall.index,
1560                               root->u.cgrp.name, 0);
1561         }
1562         root = root->next;
1563     }
1564 }
1565
1566
1567 /********************************************************************
1568  * STATIC ANALYSIS
1569  ********************************************************************/
1570
1571 /*! \brief
1572  * Marks a subtree completely dynamic or undoes such a change.
1573  *
1574  * \param     sel      Selection subtree to mark.
1575  * \param[in] bDynamic If true, the \p bStatic flag of the whole
1576  *   selection subtree is cleared. If false, the flag is restored to
1577  *   using \ref SEL_DYNAMIC.
1578  *
1579  * Does not descend into parameters of methods unless the parameters
1580  * are evaluated for each atom.
1581  */
1582 static void
1583 mark_subexpr_dynamic(const SelectionTreeElementPointer &sel,
1584                      bool bDynamic)
1585 {
1586     if (!bDynamic && !(sel->flags & SEL_DYNAMIC))
1587     {
1588         sel->cdata->flags |= SEL_CDATA_STATIC;
1589     }
1590     else
1591     {
1592         sel->cdata->flags &= ~SEL_CDATA_STATIC;
1593     }
1594     SelectionTreeElementPointer child = sel->child;
1595     while (child)
1596     {
1597         if (sel->type != SEL_EXPRESSION || child->type != SEL_SUBEXPRREF
1598             || (child->u.param->flags & SPAR_ATOMVAL))
1599         {
1600             mark_subexpr_dynamic(child, bDynamic);
1601         }
1602         child = child->next;
1603     }
1604 }
1605
1606 /*! \brief
1607  * Frees memory for subexpressions that are no longer needed.
1608  *
1609  * \param     sel      Selection subtree to check.
1610  *
1611  * Checks whether the subtree rooted at \p sel refers to any \ref SEL_SUBEXPR
1612  * elements that are not referred to by anything else except their own root
1613  * element. If such elements are found, all memory allocated for them is freed
1614  * except the actual element. The element is left because otherwise a dangling
1615  * pointer would be left at the root element, which is not traversed by this
1616  * function. Later compilation passes remove the stub elements.
1617  */
1618 static void
1619 release_subexpr_memory(const SelectionTreeElementPointer &sel)
1620 {
1621     if (sel->type == SEL_SUBEXPRREF)
1622     {
1623         const SelectionTreeElementPointer &subexpr = sel->child;
1624         if (subexpr.use_count() == 2)
1625         {
1626             release_subexpr_memory(subexpr);
1627             subexpr->name = NULL;
1628             // Free children.
1629             subexpr->child.reset();
1630             subexpr->freeValues();
1631             subexpr->freeExpressionData();
1632             subexpr->freeCompilerData();
1633         }
1634     }
1635     else
1636     {
1637         SelectionTreeElementPointer child = sel->child;
1638         while (child)
1639         {
1640             release_subexpr_memory(child);
1641             child = child->next;
1642         }
1643     }
1644 }
1645
1646 /*! \brief
1647  * Makes an evaluated selection element static.
1648  *
1649  * \param     sel   Selection element to make static.
1650  *
1651  * The evaluated value becomes the value of the static element.
1652  * The element type is changed to SEL_CONST and the children are
1653  * deleted.
1654  */
1655 static void
1656 make_static(const SelectionTreeElementPointer &sel)
1657 {
1658     /* If this is a subexpression reference and the data is stored in the
1659      * child, we transfer data ownership before doing anything else. */
1660     if (sel->type == SEL_SUBEXPRREF
1661         && (sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
1662     {
1663         if (sel->child->child->flags & SEL_ALLOCDATA)
1664         {
1665             sel->flags |= SEL_ALLOCDATA;
1666             sel->child->child->flags &= ~SEL_ALLOCDATA;
1667         }
1668         if (sel->child->child->flags & SEL_ALLOCVAL)
1669         {
1670             sel->flags |= SEL_ALLOCVAL;
1671             sel->v.nalloc = sel->child->child->v.nalloc;
1672             sel->child->child->flags &= ~SEL_ALLOCVAL;
1673             sel->child->child->v.nalloc = -1;
1674         }
1675     }
1676     /* Free the children. */
1677     release_subexpr_memory(sel);
1678     sel->child.reset();
1679     /* Free the expression data as it is no longer needed */
1680     sel->freeExpressionData();
1681     /* Make the item static */
1682     sel->name            = NULL;
1683     sel->type            = SEL_CONST;
1684     sel->evaluate        = NULL;
1685     sel->cdata->evaluate = NULL;
1686     /* Set the group value.
1687      * freeExpressionData() frees the cgrp group, so we can just override it.
1688      * */
1689     if (sel->v.type == GROUP_VALUE)
1690     {
1691         gmx_ana_index_set(&sel->u.cgrp, sel->v.u.g->isize, sel->v.u.g->index, NULL, 0);
1692     }
1693 }
1694
1695 /*! \brief
1696  * Evaluates a constant expression during analyze_static().
1697  *
1698  * \param[in]     data Evaluation data.
1699  * \param[in,out] sel Selection to process.
1700  * \param[in]     g   The evaluation group.
1701  * \returns       0 on success, a non-zero error code on error.
1702  */
1703 static void
1704 process_const(gmx_sel_evaluate_t *data,
1705               const SelectionTreeElementPointer &sel,
1706               gmx_ana_index_t *g)
1707 {
1708     if (sel->v.type == GROUP_VALUE)
1709     {
1710         if (sel->cdata->evaluate)
1711         {
1712             sel->cdata->evaluate(data, sel, g);
1713         }
1714     }
1715     /* Other constant expressions do not need evaluation */
1716 }
1717
1718 /*! \brief
1719  * Sets the parameter value pointer for \ref SEL_SUBEXPRREF params.
1720  *
1721  * \param[in,out] sel Selection to process.
1722  *
1723  * Copies the value pointer of \p sel to \c sel->u.param if one is present
1724  * and should receive the value from the compiler
1725  * (most parameter values are handled during parsing).
1726  * If \p sel is not of type \ref SEL_SUBEXPRREF, or if \c sel->u.param is NULL,
1727  * the function does nothing.
1728  * Also, if the \c sel->u.param does not have \ref SPAR_VARNUM or
1729  * \ref SPAR_ATOMVAL, the function returns immediately.
1730  */
1731 static void
1732 store_param_val(const SelectionTreeElementPointer &sel)
1733 {
1734     /* Return immediately if there is no parameter. */
1735     if (sel->type != SEL_SUBEXPRREF || !sel->u.param)
1736     {
1737         return;
1738     }
1739
1740     /* Or if the value does not need storing. */
1741     if (!(sel->u.param->flags & (SPAR_VARNUM | SPAR_ATOMVAL)))
1742     {
1743         return;
1744     }
1745
1746     if (sel->v.type == INT_VALUE || sel->v.type == REAL_VALUE
1747         || sel->v.type == STR_VALUE)
1748     {
1749         _gmx_selvalue_setstore(&sel->u.param->val, sel->v.u.ptr);
1750     }
1751 }
1752
1753 /*! \brief
1754  * Handles the initialization of a selection method during analyze_static() pass.
1755  *
1756  * \param[in,out] sel Selection element to process.
1757  * \param[in]     top Topology structure.
1758  * \param[in]     isize Size of the evaluation group for the element.
1759  * \returns       0 on success, a non-zero error code on return.
1760  *
1761  * Calls sel_initfunc() (and possibly sel_outinitfunc()) to initialize the
1762  * method.
1763  * If no \ref SPAR_ATOMVAL parameters are present, multiple initialization
1764  * is prevented by using \ref SEL_METHODINIT and \ref SEL_OUTINIT flags.
1765  */
1766 static void
1767 init_method(const SelectionTreeElementPointer &sel, t_topology *top, int isize)
1768 {
1769     /* Find out whether there are any atom-valued parameters */
1770     bool bAtomVal = false;
1771     SelectionTreeElementPointer child = sel->child;
1772     while (child)
1773     {
1774         if (child->flags & SEL_ATOMVAL)
1775         {
1776             bAtomVal = true;
1777         }
1778         child = child->next;
1779     }
1780
1781     /* Initialize the method */
1782     if (sel->u.expr.method->init
1783         && (bAtomVal || !(sel->flags & SEL_METHODINIT)))
1784     {
1785         sel->flags |= SEL_METHODINIT;
1786         sel->u.expr.method->init(top, sel->u.expr.method->nparams,
1787                 sel->u.expr.method->param, sel->u.expr.mdata);
1788     }
1789     if (bAtomVal || !(sel->flags & SEL_OUTINIT))
1790     {
1791         sel->flags |= SEL_OUTINIT;
1792         if (sel->u.expr.method->outinit)
1793         {
1794             sel->u.expr.method->outinit(top, &sel->v, sel->u.expr.mdata);
1795             if (sel->v.type != POS_VALUE && sel->v.type != GROUP_VALUE)
1796             {
1797                 alloc_selection_data(sel, isize, true);
1798             }
1799         }
1800         else
1801         {
1802             alloc_selection_data(sel, isize, true);
1803             if ((sel->flags & SEL_DYNAMIC)
1804                 && sel->v.type != GROUP_VALUE && sel->v.type != POS_VALUE)
1805             {
1806                 sel->v.nr = isize;
1807             }
1808             /* If the method is char-valued, pre-allocate the strings. */
1809             if (sel->u.expr.method->flags & SMETH_CHARVAL)
1810             {
1811                 int  i;
1812
1813                 /* A sanity check */
1814                 if (sel->v.type != STR_VALUE)
1815                 {
1816                     GMX_THROW(gmx::InternalError("Char-valued selection method in non-string element"));
1817                 }
1818                 sel->flags |= SEL_ALLOCDATA;
1819                 for (i = 0; i < isize; ++i)
1820                 {
1821                     if (sel->v.u.s[i] == NULL)
1822                     {
1823                         snew(sel->v.u.s[i], 2);
1824                     }
1825                 }
1826             }
1827         }
1828         /* Clear the values for dynamic output to avoid valgrind warnings. */
1829         if ((sel->flags & SEL_DYNAMIC) && sel->v.type == REAL_VALUE)
1830         {
1831             int i;
1832
1833             for (i = 0; i < sel->v.nr; ++i)
1834             {
1835                 sel->v.u.r[i] = 0.0;
1836             }
1837         }
1838     }
1839 }
1840
1841 /*! \brief
1842  * Evaluates the static part of a boolean expression.
1843  *
1844  * \param[in]     data Evaluation data.
1845  * \param[in,out] sel Boolean selection element whose children should be
1846  *   processed.
1847  * \param[in]     g   The evaluation group.
1848  * \returns       0 on success, a non-zero error code on error.
1849  *
1850  * reorder_item_static_children() should have been called.
1851  */
1852 static void
1853 evaluate_boolean_static_part(gmx_sel_evaluate_t *data,
1854                              const SelectionTreeElementPointer &sel,
1855                              gmx_ana_index_t *g)
1856 {
1857     /* Find the last static subexpression */
1858     SelectionTreeElementPointer child = sel->child;
1859     while (child->next && (child->next->cdata->flags & SEL_CDATA_STATIC))
1860     {
1861         child = child->next;
1862     }
1863     if (!(child->cdata->flags & SEL_CDATA_STATIC))
1864     {
1865         return;
1866     }
1867
1868     /* Evalute the static part if there is more than one expression */
1869     if (child != sel->child)
1870     {
1871         SelectionTreeElementPointer next = child->next;
1872         child->next.reset();
1873         sel->cdata->evaluate(data, sel, g);
1874         /* Replace the subexpressions with the result */
1875         child.reset(new SelectionTreeElement(SEL_CONST));
1876         child->flags      = SEL_FLAGSSET | SEL_SINGLEVAL | SEL_ALLOCVAL | SEL_ALLOCDATA;
1877         _gmx_selelem_set_vtype(child, GROUP_VALUE);
1878         child->evaluate   = NULL;
1879         _gmx_selvalue_reserve(&child->v, 1);
1880         gmx_ana_index_copy(child->v.u.g, sel->v.u.g, true);
1881         init_item_compilerdata(child);
1882         init_item_minmax_groups(child);
1883         child->cdata->flags &= ~SEL_CDATA_STATICEVAL;
1884         child->cdata->flags |= sel->cdata->flags & SEL_CDATA_STATICEVAL;
1885         child->next = next;
1886         // Frees the old static subexpressions.
1887         sel->child = child;
1888     }
1889     else if (child->evaluate)
1890     {
1891         child->evaluate(data, child, g);
1892     }
1893     /* Set the evaluation function for the constant element.
1894      * We never need to evaluate the element again during compilation,
1895      * but we may need to evaluate the static part again if the
1896      * expression is not an OR with a static evaluation group.
1897      * If we reach here with a NOT expression, the NOT expression
1898      * is also static, and will be made a constant later, so don't waste
1899      * time copying the group. */
1900     child->evaluate = NULL;
1901     if (sel->u.boolt == BOOL_NOT
1902         || ((sel->cdata->flags & SEL_CDATA_STATICEVAL)
1903             && sel->u.boolt == BOOL_OR))
1904     {
1905         child->cdata->evaluate = NULL;
1906     }
1907     else
1908     {
1909         child->cdata->evaluate = &_gmx_sel_evaluate_static;
1910         /* The cgrp has only been allocated if it originated from an
1911          * external index group. In that case, we need special handling
1912          * to preserve the name of the group and to not leak memory.
1913          * If cgrp has been set in make_static(), it is not allocated,
1914          * and hence we can overwrite it safely. */
1915         if (child->u.cgrp.nalloc_index > 0)
1916         {
1917             char *name = child->u.cgrp.name;
1918             gmx_ana_index_copy(&child->u.cgrp, child->v.u.g, false);
1919             gmx_ana_index_squeeze(&child->u.cgrp);
1920             child->u.cgrp.name = name;
1921         }
1922         else
1923         {
1924             gmx_ana_index_copy(&child->u.cgrp, child->v.u.g, true);
1925         }
1926     }
1927 }
1928
1929 /*! \brief
1930  * Evaluates the minimum and maximum groups for a boolean expression.
1931  *
1932  * \param[in]  sel  \ref SEL_BOOLEAN element currently being evaluated.
1933  * \param[in]  g    Group for which \p sel has been evaluated.
1934  * \param[out] gmin Largest subset of the possible values of \p sel.
1935  * \param[out] gmax Smallest superset of the possible values of \p sel.
1936  *
1937  * This is a helper function for analyze_static() that is called for
1938  * dynamic \ref SEL_BOOLEAN elements after they have been evaluated.
1939  * It uses the minimum and maximum groups of the children to calculate
1940  * the minimum and maximum groups for \p sel, and also updates the static
1941  * part of \p sel (which is in the first child) if the children give
1942  * cause for this.
1943  *
1944  * This function may allocate some extra memory for \p gmin and \p gmax,
1945  * but as these groups are freed at the end of analyze_static() (which is
1946  * reached shortly after this function returns), this should not be a major
1947  * problem.
1948  */
1949 static void
1950 evaluate_boolean_minmax_grps(const SelectionTreeElementPointer &sel,
1951                              gmx_ana_index_t *g,
1952                              gmx_ana_index_t *gmin, gmx_ana_index_t *gmax)
1953 {
1954     SelectionTreeElementPointer child;
1955
1956     switch (sel->u.boolt)
1957     {
1958         case BOOL_NOT:
1959             gmx_ana_index_reserve(gmin, g->isize);
1960             gmx_ana_index_reserve(gmax, g->isize);
1961             gmx_ana_index_difference(gmax, g, sel->child->cdata->gmin);
1962             gmx_ana_index_difference(gmin, g, sel->child->cdata->gmax);
1963             break;
1964
1965         case BOOL_AND:
1966             gmx_ana_index_copy(gmin, sel->child->cdata->gmin, true);
1967             gmx_ana_index_copy(gmax, sel->child->cdata->gmax, true);
1968             child = sel->child->next;
1969             while (child && gmax->isize > 0)
1970             {
1971                 gmx_ana_index_intersection(gmin, gmin, child->cdata->gmin);
1972                 gmx_ana_index_intersection(gmax, gmax, child->cdata->gmax);
1973                 child = child->next;
1974             }
1975             /* Update the static part if other expressions limit it */
1976             if ((sel->child->cdata->flags & SEL_CDATA_STATIC)
1977                 && sel->child->v.u.g->isize > gmax->isize)
1978             {
1979                 gmx_ana_index_copy(sel->child->v.u.g, gmax, false);
1980                 gmx_ana_index_squeeze(sel->child->v.u.g);
1981                 if (sel->child->u.cgrp.isize > 0)
1982                 {
1983                     gmx_ana_index_copy(&sel->child->u.cgrp, gmax, false);
1984                     gmx_ana_index_squeeze(&sel->child->u.cgrp);
1985                 }
1986             }
1987             break;
1988
1989         case BOOL_OR:
1990             /* We can assume here that the gmin of children do not overlap
1991              * because of the way _gmx_sel_evaluate_or() works. */
1992             gmx_ana_index_reserve(gmin, g->isize);
1993             gmx_ana_index_reserve(gmax, g->isize);
1994             gmx_ana_index_copy(gmin, sel->child->cdata->gmin, false);
1995             gmx_ana_index_copy(gmax, sel->child->cdata->gmax, false);
1996             child = sel->child->next;
1997             while (child && gmin->isize < g->isize)
1998             {
1999                 gmx_ana_index_merge(gmin, gmin, child->cdata->gmin);
2000                 gmx_ana_index_union(gmax, gmax, child->cdata->gmax);
2001                 child = child->next;
2002             }
2003             /* Update the static part if other expressions have static parts
2004              * that are not included. */
2005             if ((sel->child->cdata->flags & SEL_CDATA_STATIC)
2006                 && sel->child->v.u.g->isize < gmin->isize)
2007             {
2008                 GMX_RELEASE_ASSERT(sel->child->type == SEL_CONST,
2009                         "The first child should have already been evaluated "
2010                         "to a constant expression");
2011                 gmx_ana_index_reserve(sel->child->v.u.g, gmin->isize);
2012                 gmx_ana_index_copy(sel->child->v.u.g, gmin, false);
2013                 if (sel->child->u.cgrp.nalloc_index > 0)
2014                 {
2015                     /* Keep the name as in evaluate_boolean_static_part(). */
2016                     char *name = sel->child->u.cgrp.name;
2017                     gmx_ana_index_reserve(&sel->child->u.cgrp, gmin->isize);
2018                     gmx_ana_index_copy(&sel->child->u.cgrp, gmin, false);
2019                     sel->child->u.cgrp.name = name;
2020                 }
2021                 else
2022                 {
2023                     GMX_RELEASE_ASSERT(sel->child->u.cgrp.index == sel->child->v.u.g->index,
2024                             "If not allocated, the static group should equal the value");
2025                     sel->child->u.cgrp.isize = sel->child->v.u.g->isize;
2026                 }
2027             }
2028             break;
2029
2030         case BOOL_XOR: /* Should not be reached */
2031             GMX_THROW(gmx::NotImplementedError("xor expressions not implemented"));
2032             break;
2033     }
2034 }
2035
2036 /*! \brief
2037  * Evaluates the static parts of \p sel and analyzes the structure.
2038  *
2039  * \param[in]     data Evaluation data.
2040  * \param[in,out] sel  Selection currently being evaluated.
2041  * \param[in]     g    Group for which \p sel should be evaluated.
2042  * \returns       0 on success, a non-zero error code on error.
2043  *
2044  * This function is used as the replacement for the
2045  * gmx::SelectionTreeElement::evaluate function pointer.
2046  * It does the single most complex task in the compiler: after all elements
2047  * have been processed, the \p gmin and \p gmax fields of \p t_compiler_data
2048  * have been properly initialized, enough memory has been allocated for
2049  * storing the value of each expression, and the static parts of the
2050  * expressions have been evaluated.
2051  * The above is exactly true only for elements other than subexpressions:
2052  * another pass is required for subexpressions that are referred to more than
2053  * once and whose evaluation group is not known in advance.
2054  */
2055 static void
2056 analyze_static(gmx_sel_evaluate_t *data,
2057                const SelectionTreeElementPointer &sel,
2058                gmx_ana_index_t *g)
2059 {
2060     bool             bDoMinMax;
2061
2062     if (sel->type != SEL_ROOT && g)
2063     {
2064         alloc_selection_data(sel, g->isize, false);
2065     }
2066
2067     bDoMinMax = (sel->cdata->flags & SEL_CDATA_DOMINMAX);
2068     if (sel->type != SEL_SUBEXPR && bDoMinMax)
2069     {
2070         gmx_ana_index_deinit(sel->cdata->gmin);
2071         gmx_ana_index_deinit(sel->cdata->gmax);
2072     }
2073
2074     /* TODO: This switch is awfully long... */
2075     switch (sel->type)
2076     {
2077         case SEL_CONST:
2078             process_const(data, sel, g);
2079             break;
2080
2081         case SEL_EXPRESSION:
2082         case SEL_MODIFIER:
2083             GMX_ASSERT(g, "group cannot be null");
2084             _gmx_sel_evaluate_method_params(data, sel, g);
2085             init_method(sel, data->top, g->isize);
2086             if (!(sel->flags & SEL_DYNAMIC))
2087             {
2088                 sel->cdata->evaluate(data, sel, g);
2089                 if (sel->cdata->flags & SEL_CDATA_STATIC)
2090                 {
2091                     make_static(sel);
2092                 }
2093             }
2094             else
2095             {
2096                 /* Modifiers need to be evaluated even though they process
2097                  * positions to get the modified output groups from the
2098                  * maximum possible selections. */
2099                 if (sel->type == SEL_MODIFIER)
2100                 {
2101                     sel->cdata->evaluate(data, sel, g);
2102                 }
2103                 if (bDoMinMax)
2104                 {
2105                     gmx_ana_index_copy(sel->cdata->gmax, g, true);
2106                 }
2107             }
2108             break;
2109
2110         case SEL_BOOLEAN:
2111             if (!(sel->flags & SEL_DYNAMIC))
2112             {
2113                 sel->cdata->evaluate(data, sel, g);
2114                 if (sel->cdata->flags & SEL_CDATA_STATIC)
2115                 {
2116                     make_static(sel);
2117                 }
2118             }
2119             else
2120             {
2121                 /* Evalute the static part if there is more than one expression */
2122                 evaluate_boolean_static_part(data, sel, g);
2123
2124                 /* Evaluate the selection.
2125                  * If the type is boolean, we must explicitly handle the
2126                  * static part evaluated in evaluate_boolean_static_part()
2127                  * here because g may be larger. */
2128                 if (sel->u.boolt == BOOL_AND && sel->child->type == SEL_CONST)
2129                 {
2130                     sel->cdata->evaluate(data, sel, sel->child->v.u.g);
2131                 }
2132                 else
2133                 {
2134                     sel->cdata->evaluate(data, sel, g);
2135                 }
2136
2137                 /* Evaluate minimal and maximal selections */
2138                 evaluate_boolean_minmax_grps(sel, g, sel->cdata->gmin,
2139                                              sel->cdata->gmax);
2140             }
2141             break;
2142
2143         case SEL_ARITHMETIC:
2144             sel->cdata->evaluate(data, sel, g);
2145             if (!(sel->flags & SEL_DYNAMIC))
2146             {
2147                 if (sel->cdata->flags & SEL_CDATA_STATIC)
2148                 {
2149                     make_static(sel);
2150                 }
2151             }
2152             else if (bDoMinMax)
2153             {
2154                 gmx_ana_index_copy(sel->cdata->gmax, g, true);
2155             }
2156             break;
2157
2158         case SEL_ROOT:
2159             sel->cdata->evaluate(data, sel, g);
2160             break;
2161
2162         case SEL_SUBEXPR:
2163             if (sel->cdata->flags & (SEL_CDATA_SIMPLESUBEXPR | SEL_CDATA_FULLEVAL))
2164             {
2165                 sel->cdata->evaluate(data, sel, g);
2166                 _gmx_selvalue_setstore(&sel->v, sel->child->v.u.ptr);
2167             }
2168             else if (sel->u.cgrp.isize == 0)
2169             {
2170                 GMX_ASSERT(g, "group cannot be null");
2171                 gmx_ana_index_reserve(&sel->u.cgrp, g->isize);
2172                 sel->cdata->evaluate(data, sel, g);
2173                 if (bDoMinMax)
2174                 {
2175                     gmx_ana_index_copy(sel->cdata->gmin, sel->child->cdata->gmin, true);
2176                     gmx_ana_index_copy(sel->cdata->gmax, sel->child->cdata->gmax, true);
2177                 }
2178             }
2179             else
2180             {
2181                 int isize = gmx_ana_index_difference_size(g, &sel->u.cgrp);
2182                 if (isize > 0)
2183                 {
2184                     isize += sel->u.cgrp.isize;
2185                     gmx_ana_index_reserve(&sel->u.cgrp, isize);
2186                     alloc_selection_data(sel, isize, false);
2187                 }
2188                 sel->cdata->evaluate(data, sel, g);
2189                 if (isize > 0 && bDoMinMax)
2190                 {
2191                     gmx_ana_index_reserve(sel->cdata->gmin,
2192                                           sel->cdata->gmin->isize
2193                                           + sel->child->cdata->gmin->isize);
2194                     gmx_ana_index_reserve(sel->cdata->gmax,
2195                                           sel->cdata->gmax->isize
2196                                           + sel->child->cdata->gmax->isize);
2197                     gmx_ana_index_merge(sel->cdata->gmin, sel->cdata->gmin,
2198                                         sel->child->cdata->gmin);
2199                     gmx_ana_index_merge(sel->cdata->gmax, sel->cdata->gmax,
2200                                         sel->child->cdata->gmax);
2201                 }
2202             }
2203             break;
2204
2205         case SEL_SUBEXPRREF:
2206             if (!g && !(sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
2207             {
2208                 /* The subexpression should have been evaluated if g is NULL
2209                  * (i.e., this is a method parameter or a direct value of a
2210                  * selection). */
2211                 alloc_selection_data(sel, sel->child->cdata->gmax->isize, true);
2212             }
2213             sel->cdata->evaluate(data, sel, g);
2214             if ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
2215                 && (sel->child->child->flags & SEL_ALLOCVAL))
2216             {
2217                 _gmx_selvalue_setstore(&sel->v, sel->child->child->v.u.ptr);
2218             }
2219             /* Store the parameter value if required */
2220             store_param_val(sel);
2221             if (!(sel->flags & SEL_DYNAMIC))
2222             {
2223                 if (sel->cdata->flags & SEL_CDATA_STATIC)
2224                 {
2225                     make_static(sel);
2226                 }
2227             }
2228             else if (bDoMinMax)
2229             {
2230                 if ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR) || !g)
2231                 {
2232                     gmx_ana_index_copy(sel->cdata->gmin, sel->child->cdata->gmin, true);
2233                     gmx_ana_index_copy(sel->cdata->gmax, sel->child->cdata->gmax, true);
2234                 }
2235                 else
2236                 {
2237                     gmx_ana_index_reserve(sel->cdata->gmin,
2238                                           min(g->isize, sel->child->cdata->gmin->isize));
2239                     gmx_ana_index_reserve(sel->cdata->gmax,
2240                                           min(g->isize, sel->child->cdata->gmax->isize));
2241                     gmx_ana_index_intersection(sel->cdata->gmin,
2242                                                sel->child->cdata->gmin, g);
2243                     gmx_ana_index_intersection(sel->cdata->gmax,
2244                                                sel->child->cdata->gmax, g);
2245                 }
2246             }
2247             break;
2248
2249         case SEL_GROUPREF: /* Should not be reached */
2250             GMX_THROW(gmx::APIError("Unresolved group reference in compilation"));
2251     }
2252
2253     /* Update the minimal and maximal evaluation groups */
2254     if (bDoMinMax)
2255     {
2256         gmx_ana_index_squeeze(sel->cdata->gmin);
2257         gmx_ana_index_squeeze(sel->cdata->gmax);
2258         sfree(sel->cdata->gmin->name);
2259         sfree(sel->cdata->gmax->name);
2260         sel->cdata->gmin->name = NULL;
2261         sel->cdata->gmax->name = NULL;
2262     }
2263
2264     /* Replace the result of the evaluation */
2265     /* This is not necessary for subexpressions or for boolean negations
2266      * because the evaluation function already has done it properly. */
2267     if (sel->v.type == GROUP_VALUE && (sel->flags & SEL_DYNAMIC)
2268         && sel->type != SEL_SUBEXPR
2269         && !(sel->type == SEL_BOOLEAN && sel->u.boolt == BOOL_NOT))
2270     {
2271         if (sel->cdata->flags & SEL_CDATA_EVALMAX)
2272         {
2273             gmx_ana_index_copy(sel->v.u.g, sel->cdata->gmax, false);
2274         }
2275         else
2276         {
2277             gmx_ana_index_copy(sel->v.u.g, sel->cdata->gmin, false);
2278         }
2279     }
2280 }
2281
2282
2283 /********************************************************************
2284  * EVALUATION GROUP INITIALIZATION
2285  ********************************************************************/
2286
2287 /*! \brief
2288  * Initializes the evaluation group for a \ref SEL_ROOT element.
2289  *
2290  * \param     root Root element to initialize.
2291  * \param[in] gall Group of all atoms.
2292  *
2293  * Checks whether it is necessary to evaluate anything through the root
2294  * element, and either clears the evaluation function or initializes the
2295  * evaluation group.
2296  */
2297 static void
2298 init_root_item(const SelectionTreeElementPointer &root,
2299                gmx_ana_index_t *gall)
2300 {
2301     const SelectionTreeElementPointer &expr = root->child;
2302     /* Subexpressions with non-static evaluation group should not be
2303      * evaluated by the root, and neither should be single-reference
2304      * subexpressions that don't evaluate for all atoms. */
2305     if (expr->type == SEL_SUBEXPR
2306         && (!(root->child->cdata->flags & SEL_CDATA_STATICEVAL)
2307             || ((root->child->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
2308                 && !(root->child->cdata->flags & SEL_CDATA_FULLEVAL))))
2309     {
2310         root->evaluate = NULL;
2311         if (root->cdata)
2312         {
2313             root->cdata->evaluate = NULL;
2314         }
2315     }
2316
2317     /* Set the evaluation group */
2318     char *name = root->u.cgrp.name;
2319     if (root->evaluate)
2320     {
2321         /* Non-atom-valued non-group expressions don't care about the group, so
2322          * don't allocate any memory for it. */
2323         if ((expr->flags & SEL_VARNUMVAL)
2324             || ((expr->flags & SEL_SINGLEVAL) && expr->v.type != GROUP_VALUE))
2325         {
2326             gmx_ana_index_set(&root->u.cgrp, -1, NULL, NULL, 0);
2327         }
2328         else if (expr->cdata->gmax->isize == gall->isize)
2329         {
2330             /* Save some memory by only referring to the global group. */
2331             gmx_ana_index_set(&root->u.cgrp, gall->isize, gall->index, NULL, 0);
2332         }
2333         else
2334         {
2335             gmx_ana_index_copy(&root->u.cgrp, expr->cdata->gmax, true);
2336         }
2337         /* For selections, store the maximum group for
2338          * gmx_ana_selcollection_evaluate_fin() as the value of the root
2339          * element (unused otherwise). */
2340         if (expr->type != SEL_SUBEXPR && expr->v.u.p->g)
2341         {
2342             SelectionTreeElementPointer child = expr;
2343
2344             /* TODO: This code is copied from parsetree.c; it would be better
2345              * to have this hardcoded only in one place. */
2346             while (child->type == SEL_MODIFIER)
2347             {
2348                 child = child->child;
2349                 if (child->type == SEL_SUBEXPRREF)
2350                 {
2351                     child = child->child->child;
2352                 }
2353             }
2354             if (child->type == SEL_SUBEXPRREF)
2355             {
2356                 child = child->child->child;
2357             }
2358             if (child->child->flags & SEL_DYNAMIC)
2359             {
2360                 _gmx_selelem_set_vtype(root, GROUP_VALUE);
2361                 root->flags  |= (SEL_ALLOCVAL | SEL_ALLOCDATA);
2362                 _gmx_selvalue_reserve(&root->v, 1);
2363                 gmx_ana_index_copy(root->v.u.g, expr->v.u.p->g, true);
2364             }
2365         }
2366     }
2367     else
2368     {
2369         gmx_ana_index_clear(&root->u.cgrp);
2370     }
2371     root->u.cgrp.name = name;
2372 }
2373
2374
2375 /********************************************************************
2376  * FINAL SUBEXPRESSION OPTIMIZATION
2377  ********************************************************************/
2378
2379 /*! \brief
2380  * Optimizes subexpression evaluation.
2381  *
2382  * \param     sel Root of the selection subtree to process.
2383  *
2384  * Optimizes away some unnecessary evaluation of subexpressions that are only
2385  * referenced once.
2386  */
2387 static void
2388 postprocess_item_subexpressions(const SelectionTreeElementPointer &sel)
2389 {
2390     GMX_ASSERT(!(sel->child == NULL &&
2391                  (sel->type == SEL_SUBEXPRREF || sel->type == SEL_SUBEXPR)),
2392                "Subexpression elements should always have a child element");
2393
2394     /* Process children. */
2395     if (sel->type != SEL_SUBEXPRREF)
2396     {
2397         SelectionTreeElementPointer child = sel->child;
2398         while (child)
2399         {
2400             postprocess_item_subexpressions(child);
2401             child = child->next;
2402         }
2403     }
2404
2405     /* Replace the evaluation function of statically evaluated subexpressions
2406      * for which the static group was not known in advance. */
2407     if (sel->type == SEL_SUBEXPR && sel->cdata->refcount > 1
2408         && (sel->cdata->flags & SEL_CDATA_STATICEVAL)
2409         && !(sel->cdata->flags & SEL_CDATA_FULLEVAL))
2410     {
2411         char *name;
2412
2413         /* We need to free memory allocated for the group, because it is no
2414          * longer needed (and would be lost on next call to the evaluation
2415          * function). But we need to preserve the name. */
2416         name = sel->u.cgrp.name;
2417         gmx_ana_index_deinit(&sel->u.cgrp);
2418         sel->u.cgrp.name = name;
2419
2420         sel->evaluate = &_gmx_sel_evaluate_subexpr_staticeval;
2421         sel->cdata->evaluate = sel->evaluate;
2422
2423         sel->child->freeValues();
2424         sel->child->mempool = NULL;
2425         _gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
2426         sel->child->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
2427     }
2428
2429     /* Adjust memory allocation flags for subexpressions that are used only
2430      * once.  This is not strictly necessary, but we do it to have the memory
2431      * managed consistently for all types of subexpressions. */
2432     if (sel->type == SEL_SUBEXPRREF
2433         && (sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR))
2434     {
2435         if (sel->child->child->flags & SEL_ALLOCVAL)
2436         {
2437             sel->flags |= SEL_ALLOCVAL;
2438             sel->flags |= (sel->child->child->flags & SEL_ALLOCDATA);
2439             sel->v.nalloc = sel->child->child->v.nalloc;
2440             sel->child->child->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
2441             sel->child->child->v.nalloc = -1;
2442         }
2443     }
2444
2445     /* Do the same for subexpressions that are evaluated at once for all atoms. */
2446     if (sel->type == SEL_SUBEXPR
2447         && !(sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
2448         && (sel->cdata->flags & SEL_CDATA_FULLEVAL))
2449     {
2450         sel->flags |= SEL_ALLOCVAL;
2451         sel->flags |= (sel->child->flags & SEL_ALLOCDATA);
2452         sel->v.nalloc = sel->child->v.nalloc;
2453         sel->child->flags &= ~(SEL_ALLOCVAL | SEL_ALLOCDATA);
2454         sel->child->v.nalloc = -1;
2455     }
2456 }
2457
2458
2459 /********************************************************************
2460  * COM CALCULATION INITIALIZATION
2461  ********************************************************************/
2462
2463 /*! \brief
2464  * Initializes COM/COG calculation for method expressions that require it.
2465  *
2466  * \param     sel    Selection subtree to process.
2467  * \param[in,out] pcc   Position calculation collection to use.
2468  * \param[in] type   Default position calculation type.
2469  * \param[in] flags  Flags for default position calculation.
2470  *
2471  * Searches recursively through the selection tree for dynamic
2472  * \ref SEL_EXPRESSION elements that define the \c gmx_ana_selmethod_t::pupdate
2473  * function.
2474  * For each such element found, position calculation is initialized
2475  * for the maximal evaluation group.
2476  * The type of the calculation is determined by \p type and \p flags.
2477  * No calculation is initialized if \p type equals \ref POS_ATOM and
2478  * the method also defines the \c gmx_ana_selmethod_t::update method.
2479  */
2480 static void
2481 init_item_comg(const SelectionTreeElementPointer &sel,
2482                gmx::PositionCalculationCollection *pcc,
2483                e_poscalc_t type, int flags)
2484 {
2485     /* Initialize COM calculation for dynamic selections now that we know the maximal evaluation group */
2486     if (sel->type == SEL_EXPRESSION && sel->u.expr.method
2487         && sel->u.expr.method->pupdate)
2488     {
2489         if (!sel->u.expr.method->update || type != POS_ATOM)
2490         {
2491             /* Create a default calculation if one does not yet exist */
2492             int cflags = 0;
2493             if (!(sel->cdata->flags & SEL_CDATA_STATICEVAL))
2494             {
2495                 cflags |= POS_DYNAMIC;
2496             }
2497             if (!sel->u.expr.pc)
2498             {
2499                 cflags |= flags;
2500                 sel->u.expr.pc = pcc->createCalculation(type, cflags);
2501             }
2502             else
2503             {
2504                 gmx_ana_poscalc_set_flags(sel->u.expr.pc, cflags);
2505             }
2506             gmx_ana_poscalc_set_maxindex(sel->u.expr.pc, sel->cdata->gmax);
2507             snew(sel->u.expr.pos, 1);
2508             gmx_ana_poscalc_init_pos(sel->u.expr.pc, sel->u.expr.pos);
2509         }
2510     }
2511
2512     /* Call recursively for all children unless the children have already been processed */
2513     if (sel->type != SEL_SUBEXPRREF)
2514     {
2515         SelectionTreeElementPointer child = sel->child;
2516         while (child)
2517         {
2518             init_item_comg(child, pcc, type, flags);
2519             child = child->next;
2520         }
2521     }
2522 }
2523
2524
2525 /********************************************************************
2526  * COMPILER DATA FREEING
2527  ********************************************************************/
2528
2529 /*! \brief
2530  * Frees the allocated compiler data recursively.
2531  *
2532  * \param     sel Root of the selection subtree to process.
2533  *
2534  * Frees the data allocated for the compilation process.
2535  */
2536 static void
2537 free_item_compilerdata(const SelectionTreeElementPointer &sel)
2538 {
2539     /* Free compilation data */
2540     sel->freeCompilerData();
2541
2542     /* Call recursively for all children unless the children have already been processed */
2543     if (sel->type != SEL_SUBEXPRREF)
2544     {
2545         SelectionTreeElementPointer child = sel->child;
2546         while (child)
2547         {
2548             free_item_compilerdata(child);
2549             child = child->next;
2550         }
2551     }
2552 }
2553
2554
2555 /********************************************************************
2556  * MAIN COMPILATION FUNCTION
2557  ********************************************************************/
2558
2559 namespace gmx
2560 {
2561
2562 SelectionCompiler::SelectionCompiler()
2563 {
2564 }
2565
2566 /*!
2567  * \param[in,out] coll Selection collection to be compiled.
2568  * \returns       0 on successful compilation, a non-zero error code on error.
2569  *
2570  * Before compilation, the selection collection should have been initialized
2571  * with gmx_ana_selcollection_parse_*().
2572  * The compiled selection collection can be passed to
2573  * gmx_ana_selcollection_evaluate() to evaluate the selection for a frame.
2574  * If an error occurs, \p sc is cleared.
2575  *
2576  * The covered fraction information in \p sc is initialized to
2577  * \ref CFRAC_NONE.
2578  */
2579 void
2580 SelectionCompiler::compile(SelectionCollection *coll)
2581 {
2582     gmx_ana_selcollection_t *sc = &coll->impl_->sc_;
2583     gmx_sel_evaluate_t  evaldata;
2584     SelectionTreeElementPointer item;
2585     e_poscalc_t  post;
2586     size_t       i;
2587     int          flags;
2588     bool         bDebug = (coll->impl_->debugLevel_ >= 2
2589                            && coll->impl_->debugLevel_ != 3);
2590
2591     /* FIXME: Clean up the collection on exceptions */
2592
2593     sc->mempool = _gmx_sel_mempool_create();
2594     _gmx_sel_evaluate_init(&evaldata, sc->mempool, &sc->gall,
2595                            sc->top, NULL, NULL);
2596
2597     /* Clear the symbol table because it is not possible to parse anything
2598      * after compilation, and variable references in the symbol table can
2599      * also mess up the compilation and/or become invalid.
2600      */
2601     coll->impl_->clearSymbolTable();
2602
2603     /* Loop through selections and initialize position keyword defaults if no
2604      * other value has been provided.
2605      */
2606     for (i = 0; i < sc->sel.size(); ++i)
2607     {
2608         gmx::internal::SelectionData &sel = *sc->sel[i];
2609         init_pos_keyword_defaults(&sel.rootElement(),
2610                                   coll->impl_->spost_.c_str(),
2611                                   coll->impl_->rpost_.c_str(),
2612                                   &sel);
2613     }
2614
2615     /* Remove any unused variables. */
2616     sc->root = remove_unused_subexpressions(sc->root);
2617     /* Extract subexpressions into separate roots */
2618     sc->root = extract_subexpressions(sc->root);
2619
2620     /* Initialize the evaluation callbacks and process the tree structure
2621      * to conform to the expectations of the callback functions. */
2622     /* Also, initialize and allocate the compiler data structure */
2623     item = sc->root;
2624     while (item)
2625     {
2626         /* Process boolean and arithmetic expressions. */
2627         optimize_boolean_expressions(item);
2628         reorder_boolean_static_children(item);
2629         optimize_arithmetic_expressions(item);
2630         /* Initialize the compiler data */
2631         init_item_compilerdata(item);
2632         init_item_staticeval(item);
2633         init_item_subexpr_refcount(item);
2634         item = item->next;
2635     }
2636     /* Initialize subexpression flags.
2637      * Requires compiler flags for the full tree. */
2638     item = sc->root;
2639     while (item)
2640     {
2641         init_item_subexpr_flags(item);
2642         item = item->next;
2643     }
2644     /* Initialize evaluation.
2645      * Requires subexpression flags. */
2646     item = sc->root;
2647     while (item)
2648     {
2649         init_item_evalfunc(item);
2650         setup_memory_pooling(item, sc->mempool);
2651         init_item_evaloutput(item);
2652         item = item->next;
2653     }
2654     /* Initialize minimum/maximum index groups.
2655      * Requires evaluation output for the full tree. */
2656     item = sc->root;
2657     while (item)
2658     {
2659         init_item_minmax_groups(item);
2660         item = item->next;
2661     }
2662     /* Initialize the evaluation index groups */
2663     initialize_evalgrps(sc);
2664
2665     if (bDebug)
2666     {
2667         fprintf(stderr, "\nTree after initial compiler processing:\n");
2668         coll->printTree(stderr, false);
2669     }
2670
2671     /* Evaluate all static parts of the selection and analyze the tree
2672      * to allocate enough memory to store the value of each dynamic subtree. */
2673     item = sc->root;
2674     while (item)
2675     {
2676         if (item->child->cdata->flags & SEL_CDATA_COMMONSUBEXPR)
2677         {
2678             mark_subexpr_dynamic(item->child, true);
2679         }
2680         set_evaluation_function(item, &analyze_static);
2681         item->evaluate(&evaldata, item, NULL);
2682         item = item->next;
2683     }
2684
2685     /* At this point, static subexpressions no longer have references to them,
2686      * so they can be removed. */
2687     sc->root = remove_unused_subexpressions(sc->root);
2688     // Update the reference counts for consistency (only used for the
2689     // debugging output below).
2690     item = sc->root;
2691     while (item)
2692     {
2693         init_item_subexpr_refcount(item);
2694         item = item->next;
2695     }
2696
2697     if (bDebug)
2698     {
2699         fprintf(stderr, "\nTree after first analysis pass:\n");
2700         coll->printTree(stderr, false);
2701     }
2702
2703     /* Do a second pass to evaluate static parts of common subexpressions */
2704     item = sc->root;
2705     while (item)
2706     {
2707         if (item->child->cdata->flags & SEL_CDATA_COMMONSUBEXPR)
2708         {
2709             bool bMinMax = item->child->cdata->flags & SEL_CDATA_DOMINMAX;
2710
2711             mark_subexpr_dynamic(item->child, false);
2712             item->child->u.cgrp.isize = 0;
2713             /* We won't clear item->child->v.u.g here, because it may
2714              * be static, and hence actually point to item->child->cdata->gmax,
2715              * which is used below. We could also check whether this is the
2716              * case and only clear the group otherwise, but because the value
2717              * is actually overwritten immediately in the evaluate call, we
2718              * won't, because similar problems may arise if gmax handling ever
2719              * changes and the check were not updated.
2720              * For the same reason, we clear the min/max flag so that the
2721              * evaluation group doesn't get messed up. */
2722             set_evaluation_function(item, &analyze_static);
2723             item->child->cdata->flags &= ~SEL_CDATA_DOMINMAX;
2724             item->evaluate(&evaldata, item->child, item->child->cdata->gmax);
2725             if (bMinMax)
2726             {
2727                 item->child->cdata->flags |= SEL_CDATA_DOMINMAX;
2728             }
2729         }
2730         item = item->next;
2731     }
2732
2733     /* We need a yet another pass of subexpression removal to remove static
2734      * subexpressions referred to by common dynamic subexpressions. */
2735     sc->root = remove_unused_subexpressions(sc->root);
2736     // Update the reference counts, used by postprocess_item_subexpressions().
2737     item = sc->root;
2738     while (item)
2739     {
2740         init_item_subexpr_refcount(item);
2741         item = item->next;
2742     }
2743
2744     if (bDebug)
2745     {
2746         fprintf(stderr, "\nTree after second analysis pass:\n");
2747         coll->printTree(stderr, false);
2748     }
2749
2750     /* Initialize evaluation groups, position calculations for methods, perform
2751      * some final optimization, and free the memory allocated for the
2752      * compilation. */
2753     /* By default, use whole residues/molecules. */
2754     flags = POS_COMPLWHOLE;
2755     PositionCalculationCollection::typeFromEnum(coll->impl_->rpost_.c_str(),
2756                                                 &post, &flags);
2757     item = sc->root;
2758     while (item)
2759     {
2760         init_root_item(item, &sc->gall);
2761         postprocess_item_subexpressions(item);
2762         init_item_comg(item, &sc->pcc, post, flags);
2763         free_item_compilerdata(item);
2764         item = item->next;
2765     }
2766
2767     /* Allocate memory for the evaluation memory pool. */
2768     _gmx_sel_mempool_reserve(sc->mempool, 0);
2769
2770     /* Finish up by calculating total masses and charges. */
2771     for (i = 0; i < sc->sel.size(); ++i)
2772     {
2773         sc->sel[i]->initializeMassesAndCharges(sc->top);
2774     }
2775 }
2776
2777 } // namespace gmx