src/gmxlib/selection/compiler.c

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