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