src/gmxlib/selection/compiler.c

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