src/gromacs/selection/compiler.cpp

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