src/gromacs/selection/compiler.cpp

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