src/gromacs/selection/compiler.cpp

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