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