src/gromacs/selection/compiler.cpp

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