src/gromacs/selection/compiler.cpp

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