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