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