src/gromacs/selection/evaluate.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements functions in evaluate.h.
  34  *
  35  * \todo
  36  * One of the major bottlenecks for selection performance is that all the
  37  * evaluation is carried out for atoms.
  38  * There are several cases when the evaluation could be done for residues
  39  * or molecules instead, including keywords that select by residue and
  40  * cases where residue centers are used as reference positions.
  41  * Implementing this would require a mechanism for recognizing whether
  42  * something can be evaluated by residue/molecule instead by atom, and
  43  * converting selections by residue/molecule into selections by atom
  44  * when necessary.
  45  *
  46  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  47  * \ingroup module_selection
  48  */
  49 #include <string.h>
  50
  51 #include "gromacs/legacyheaders/maths.h"
  52 #include "gromacs/legacyheaders/smalloc.h"
  53 #include "gromacs/legacyheaders/vec.h"
  54
  55 #include "gromacs/selection/indexutil.h"
  56 #include "gromacs/selection/poscalc.h"
  57 #include "gromacs/selection/selection.h"
  58 #include "gromacs/selection/selmethod.h"
  59 #include "gromacs/utility/exceptions.h"
  60 #include "gromacs/utility/gmxassert.h"
  61
  62 #include "evaluate.h"
  63 #include "mempool.h"
  64 #include "selectioncollection-impl.h"
  65 #include "selelem.h"
  66
  67 using gmx::SelectionTreeElement;
  68 using gmx::SelectionTreeElementPointer;
  69
  70 namespace
  71 {
  72
  73 /*! \internal \brief
  74  * Reserves memory for a selection element from the evaluation memory pool.
  75  *
  76  * This class implements RAII semantics for allocating memory for selection
  77  * element values from a selection evaluation memory pool.
  78  *
  79  * \ingroup module_selection
  80  */
  81 class MempoolSelelemReserver
  82 {
  83     public:
  84         //! Constructs a reserver without initial reservation.
  85         MempoolSelelemReserver() {}
  86         /*! \brief
  87          * Constructs a reserver with initial reservation.
  88          *
  89          * \param[in,out] sel    Selection element for which to reserve.
  90          * \param[in]     count  Number of values to reserve.
  91          *
  92          * \see reserve()
  93          */
  94         MempoolSelelemReserver(const SelectionTreeElementPointer &sel, int count)
  95         {
  96             reserve(sel, count);
  97         }
  98         //! Frees any memory allocated using this reserver.
  99         ~MempoolSelelemReserver()
 100         {
 101             if (sel_)
 102             {
 103                 sel_->mempoolRelease();
 104             }
 105         }
 106
 107         /*! \brief
 108          * Reserves memory for selection element values using this reserver.
 109          *
 110          * \param[in,out] sel    Selection element for which to reserve.
 111          * \param[in]     count  Number of values to reserve.
 112          *
 113          * Allocates space to store \p count output values in \p sel from the
 114          * memory pool associated with \p sel, or from the heap if there is no
 115          * memory pool.  Type of values to allocate is automatically determined
 116          * from \p sel.
 117          */
 118         void reserve(const SelectionTreeElementPointer &sel, int count)
 119         {
 120             GMX_RELEASE_ASSERT(!sel_,
 121                                "Can only reserve one element with one instance");
 122             sel->mempoolReserve(count);
 123             sel_ = sel;
 124         }
 125
 126     private:
 127         SelectionTreeElementPointer     sel_;
 128 };
 129
 130 /*! \internal \brief
 131  * Reserves memory for an index group from the evaluation memory pool.
 132  *
 133  * This class implements RAII semantics for allocating memory for an index
 134  * group from a selection evaluation memory pool.
 135  *
 136  * \ingroup module_selection
 137  */
 138 class MempoolGroupReserver
 139 {
 140     public:
 141         /*! \brief
 142          * Creates a reserver associated with a given memory pool.
 143          *
 144          * \param    mp  Memory pool from which to reserve memory.
 145          */
 146         explicit MempoolGroupReserver(gmx_sel_mempool_t *mp)
 147             : mp_(mp), g_(NULL)
 148         {
 149         }
 150         //! Frees any memory allocated using this reserver.
 151         ~MempoolGroupReserver()
 152         {
 153             if (g_ != NULL)
 154             {
 155                 _gmx_sel_mempool_free_group(mp_, g_);
 156             }
 157         }
 158
 159         /*! \brief
 160          * Reserves memory for an index group using this reserver.
 161          *
 162          * \param[in,out] g      Index group to reserve.
 163          * \param[in]     count  Number of atoms to reserve space for.
 164          *
 165          * Allocates memory from the memory pool to store \p count atoms in
 166          * \p g.
 167          */
 168         void reserve(gmx_ana_index_t *g, int count)
 169         {
 170             GMX_RELEASE_ASSERT(g_ == NULL, "Can only reserve one element with one instance");
 171             _gmx_sel_mempool_alloc_group(mp_, g, count);
 172             g_ = g;
 173         }
 174
 175     private:
 176         gmx_sel_mempool_t      *mp_;
 177         gmx_ana_index_t        *g_;
 178 };
 179
 180 /*! \internal \brief
 181  * Assigns a temporary value for a selection element.
 182  *
 183  * This class implements RAII semantics for temporarily assigning the value
 184  * pointer of a selection element to point to a different location.
 185  *
 186  * \ingroup module_selection
 187  */
 188 class SelelemTemporaryValueAssigner
 189 {
 190     public:
 191         //! Constructs an assigner without an initial assignment.
 192         SelelemTemporaryValueAssigner()
 193             : old_ptr_(NULL), old_nalloc_(0)
 194         {
 195         }
 196         /*! \brief
 197          * Constructs an assigner with an initial assignment.
 198          *
 199          * \param[in,out] sel     Selection element for which to assign.
 200          * \param[in]     vsource Element to which \p sel values will point to.
 201          *
 202          * \see assign()
 203          */
 204         SelelemTemporaryValueAssigner(const SelectionTreeElementPointer &sel,
 205                                       const SelectionTreeElement        &vsource)
 206         {
 207             assign(sel, vsource);
 208         }
 209         //! Undoes any temporary assignment done using this assigner.
 210         ~SelelemTemporaryValueAssigner()
 211         {
 212             if (sel_)
 213             {
 214                 _gmx_selvalue_setstore_alloc(&sel_->v, old_ptr_, old_nalloc_);
 215             }
 216         }
 217
 218         /*! \brief
 219          * Assigns a temporary value pointer.
 220          *
 221          * \param[in,out] sel     Selection element for which to assign.
 222          * \param[in]     vsource Element to which \p sel values will point to.
 223          *
 224          * Assigns the value pointer in \p sel to point to the values in
 225          * \p vsource, i.e., any access/modification to values in \p sel
 226          * actually accesses values in \p vsource.
 227          */
 228         void assign(const SelectionTreeElementPointer &sel,
 229                     const SelectionTreeElement        &vsource)
 230         {
 231             GMX_RELEASE_ASSERT(!sel_,
 232                                "Can only assign one element with one instance");
 233             GMX_RELEASE_ASSERT(sel->v.type == vsource.v.type,
 234                                "Mismatching selection value types");
 235             old_ptr_    = sel->v.u.ptr;
 236             old_nalloc_ = sel->v.nalloc;
 237             _gmx_selvalue_setstore(&sel->v, vsource.v.u.ptr);
 238             sel_ = sel;
 239         }
 240
 241     private:
 242         SelectionTreeElementPointer     sel_;
 243         void                           *old_ptr_;
 244         int                             old_nalloc_;
 245 };
 246
 247 } // namespace
 248
 249 /*!
 250  * \param[in] fp       File handle to receive the output.
 251  * \param[in] evalfunc Function pointer to print.
 252  */
 253 void
 254 _gmx_sel_print_evalfunc_name(FILE *fp, gmx::sel_evalfunc evalfunc)
 255 {
 256     if (!evalfunc)
 257     {
 258         fprintf(fp, "none");
 259     }
 260     else if (evalfunc == &_gmx_sel_evaluate_root)
 261     {
 262         fprintf(fp, "root");
 263     }
 264     else if (evalfunc == &_gmx_sel_evaluate_static)
 265     {
 266         fprintf(fp, "static");
 267     }
 268     else if (evalfunc == &_gmx_sel_evaluate_subexpr_simple)
 269     {
 270         fprintf(fp, "subexpr_simple");
 271     }
 272     else if (evalfunc == &_gmx_sel_evaluate_subexpr_staticeval)
 273     {
 274         fprintf(fp, "subexpr_staticeval");
 275     }
 276     else if (evalfunc == &_gmx_sel_evaluate_subexpr)
 277     {
 278         fprintf(fp, "subexpr");
 279     }
 280     else if (evalfunc == &_gmx_sel_evaluate_subexprref_simple)
 281     {
 282         fprintf(fp, "ref_simple");
 283     }
 284     else if (evalfunc == &_gmx_sel_evaluate_subexprref)
 285     {
 286         fprintf(fp, "ref");
 287     }
 288     else if (evalfunc == &_gmx_sel_evaluate_method)
 289     {
 290         fprintf(fp, "method");
 291     }
 292     else if (evalfunc == &_gmx_sel_evaluate_modifier)
 293     {
 294         fprintf(fp, "mod");
 295     }
 296     else if (evalfunc == &_gmx_sel_evaluate_not)
 297     {
 298         fprintf(fp, "not");
 299     }
 300     else if (evalfunc == &_gmx_sel_evaluate_and)
 301     {
 302         fprintf(fp, "and");
 303     }
 304     else if (evalfunc == &_gmx_sel_evaluate_or)
 305     {
 306         fprintf(fp, "or");
 307     }
 308     else if (evalfunc == &_gmx_sel_evaluate_arithmetic)
 309     {
 310         fprintf(fp, "arithmetic");
 311     }
 312     else
 313     {
 314         fprintf(fp, "%p", (void*)(evalfunc));
 315     }
 316 }
 317
 318 /*!
 319  * \param[out] data Evaluation data structure to initialize.
 320  * \param[in]  mp   Memory pool for intermediate evaluation values.
 321  * \param[in]  gall Index group with all the atoms.
 322  * \param[in]  top  Topology structure for evaluation.
 323  * \param[in]  fr   New frame for evaluation.
 324  * \param[in]  pbc  New PBC information for evaluation.
 325  */
 326 void
 327 _gmx_sel_evaluate_init(gmx_sel_evaluate_t *data,
 328                        gmx_sel_mempool_t *mp, gmx_ana_index_t *gall,
 329                        t_topology *top, t_trxframe *fr, t_pbc *pbc)
 330 {
 331     data->mp   = mp;
 332     data->gall = gall;
 333     data->top  = top;
 334     data->fr   = fr;
 335     data->pbc  = pbc;
 336 }
 337
 338 /*! \brief
 339  * Recursively initializes the flags for evaluation.
 340  *
 341  * \param[in,out] sel Selection element to clear.
 342  *
 343  * The \ref SEL_INITFRAME flag is set for \ref SEL_EXPRESSION elements whose
 344  * method defines the \p init_frame callback (see sel_framefunc()), and
 345  * cleared for other elements.
 346  *
 347  * The \ref SEL_EVALFRAME flag is cleared for all elements.
 348  */
 349 static void
 350 init_frame_eval(SelectionTreeElementPointer sel)
 351 {
 352     while (sel)
 353     {
 354         sel->flags &= ~(SEL_INITFRAME | SEL_EVALFRAME);
 355         if (sel->type == SEL_EXPRESSION)
 356         {
 357             if (sel->u.expr.method && sel->u.expr.method->init_frame)
 358             {
 359                 sel->flags |= SEL_INITFRAME;
 360             }
 361         }
 362         if (sel->child && sel->type != SEL_SUBEXPRREF)
 363         {
 364             init_frame_eval(sel->child);
 365         }
 366         sel = sel->next;
 367     }
 368 }
 369
 370 namespace gmx
 371 {
 372
 373 SelectionEvaluator::SelectionEvaluator()
 374 {
 375 }
 376
 377 /*!
 378  * \param[in,out] coll  The selection collection to evaluate.
 379  * \param[in] fr  Frame for which the evaluation should be carried out.
 380  * \param[in] pbc PBC data, or NULL if no PBC should be used.
 381  * \returns   0 on successful evaluation, a non-zero error code on error.
 382  *
 383  * This functions sets the global variables for topology, frame and PBC,
 384  * clears some information in the selection to initialize the evaluation
 385  * for a new frame, and evaluates \p sel and all the selections pointed by
 386  * the \p next pointers of \p sel.
 387  *
 388  * This is the only function that user code should call if they want to
 389  * evaluate a selection for a new frame.
 390  */
 391 void
 392 SelectionEvaluator::evaluate(SelectionCollection *coll,
 393                              t_trxframe *fr, t_pbc *pbc)
 394 {
 395     gmx_ana_selcollection_t *sc = &coll->impl_->sc_;
 396     gmx_sel_evaluate_t       data;
 397
 398     _gmx_sel_evaluate_init(&data, sc->mempool, &sc->gall, sc->top, fr, pbc);
 399     init_frame_eval(sc->root);
 400     SelectionTreeElementPointer sel = sc->root;
 401     while (sel)
 402     {
 403         /* Clear the evaluation group of subexpressions */
 404         if (sel->child && sel->child->type == SEL_SUBEXPR)
 405         {
 406             sel->child->u.cgrp.isize = 0;
 407             /* Not strictly necessary, because the value will be overwritten
 408              * during first evaluation of the subexpression anyways, but we
 409              * clear the group for clarity. Note that this is _not_ done during
 410              * compilation because of some additional complexities involved
 411              * (see compiler.c), so it should not be relied upon in
 412              * _gmx_sel_evaluate_subexpr(). */
 413             if (sel->child->v.type == GROUP_VALUE)
 414             {
 415                 sel->child->v.u.g->isize = 0;
 416             }
 417         }
 418         if (sel->evaluate)
 419         {
 420             sel->evaluate(&data, sel, NULL);
 421         }
 422         sel = sel->next;
 423     }
 424     /* Update selection information */
 425     SelectionDataList::const_iterator isel;
 426     for (isel = sc->sel.begin(); isel != sc->sel.end(); ++isel)
 427     {
 428         internal::SelectionData &sel = **isel;
 429         sel.refreshMassesAndCharges(sc->top);
 430         sel.updateCoveredFractionForFrame();
 431     }
 432 }
 433
 434 /*!
 435  * \param[in,out] coll  The selection collection to evaluate.
 436  * \param[in]     nframes Total number of frames.
 437  */
 438 void
 439 SelectionEvaluator::evaluateFinal(SelectionCollection *coll, int nframes)
 440 {
 441     gmx_ana_selcollection_t          *sc = &coll->impl_->sc_;
 442
 443     SelectionDataList::const_iterator isel;
 444     for (isel = sc->sel.begin(); isel != sc->sel.end(); ++isel)
 445     {
 446         internal::SelectionData &sel = **isel;
 447         sel.restoreOriginalPositions(sc->top);
 448         sel.computeAverageCoveredFraction(nframes);
 449     }
 450 }
 451
 452 } // namespace gmx
 453
 454 /*!
 455  * \param[in] data Data for the current frame.
 456  * \param[in] sel  Selection element being evaluated.
 457  * \param[in] g    Group for which \p sel should be evaluated.
 458  * \returns   0 on success, a non-zero error code on error.
 459  *
 460  * Evaluates each child of \p sel in \p g.
 461  */
 462 void
 463 _gmx_sel_evaluate_children(gmx_sel_evaluate_t                *data,
 464                            const SelectionTreeElementPointer &sel,
 465                            gmx_ana_index_t                   *g)
 466 {
 467     SelectionTreeElementPointer child = sel->child;
 468     while (child)
 469     {
 470         if (child->evaluate)
 471         {
 472             child->evaluate(data, child, g);
 473         }
 474         child = child->next;
 475     }
 476 }
 477
 478 /*!
 479  * \param[in] data Data for the current frame.
 480  * \param[in] sel Selection element being evaluated.
 481  * \param[in] g   Group for which \p sel should be evaluated
 482  *   (not used, can be NULL).
 483  * \returns   0 on success, a non-zero error code on error.
 484  *
 485  * Evaluates the first child element in the group defined by \p sel->u.cgrp.
 486  * If \p sel->u.cgrp is empty, nothing is done.
 487  * The value of \p sel is not touched (root elements do not evaluate to
 488  * values).
 489  *
 490  * This function can be used as gmx::SelectionTreeElement::evaluate for
 491  * \ref SEL_ROOT elements.
 492  */
 493 void
 494 _gmx_sel_evaluate_root(gmx_sel_evaluate_t                *data,
 495                        const SelectionTreeElementPointer &sel,
 496                        gmx_ana_index_t                   *g)
 497 {
 498     if (sel->u.cgrp.isize == 0 || !sel->child->evaluate)
 499     {
 500         return;
 501     }
 502
 503     sel->child->evaluate(data, sel->child,
 504                          sel->u.cgrp.isize < 0 ? NULL : &sel->u.cgrp);
 505 }
 506
 507 /*!
 508  * \param[in] data Data for the current frame.
 509  * \param[in] sel Selection element being evaluated.
 510  * \param[in] g   Group for which \p sel should be evaluated.
 511  * \returns   0 for success.
 512  *
 513  * Sets the value of \p sel to the intersection of \p g and \p sel->u.cgrp.
 514  *
 515  * This function can be used as gmx::SelectionTreeElement::evaluate for
 516  * \ref SEL_CONST elements with value type \ref GROUP_VALUE.
 517  */
 518 void
 519 _gmx_sel_evaluate_static(gmx_sel_evaluate_t                *data,
 520                          const SelectionTreeElementPointer &sel,
 521                          gmx_ana_index_t                   *g)
 522 {
 523     gmx_ana_index_intersection(sel->v.u.g, &sel->u.cgrp, g);
 524 }
 525
 526
 527 /*********************************************************************
 528  * SUBEXPRESSION EVALUATION
 529  *********************************************************************/
 530
 531 /*!
 532  * \param[in] data Data for the current frame.
 533  * \param[in] sel  Selection element being evaluated.
 534  * \param[in] g    Group for which \p sel should be evaluated.
 535  * \returns   0 on success, a non-zero error code on error.
 536  *
 537  * Evaluates the child element (there should be exactly one) in \p g.
 538  * The compiler has taken care that the child actually stores the evaluated
 539  * value in the value pointer of this element.
 540  *
 541  * This function is used as gmx::SelectionTreeElement::evaluate for
 542  * \ref SEL_SUBEXPR elements that are used only once, and hence do not need
 543  * full subexpression handling.
 544  */
 545 void
 546 _gmx_sel_evaluate_subexpr_simple(gmx_sel_evaluate_t                *data,
 547                                  const SelectionTreeElementPointer &sel,
 548                                  gmx_ana_index_t                   *g)
 549 {
 550     if (sel->child->evaluate)
 551     {
 552         sel->child->evaluate(data, sel->child, g);
 553     }
 554     sel->v.nr = sel->child->v.nr;
 555 }
 556
 557 /*!
 558  * \param[in] data Data for the current frame.
 559  * \param[in] sel  Selection element being evaluated.
 560  * \param[in] g    Group for which \p sel should be evaluated.
 561  * \returns   0 on success, a non-zero error code on error.
 562  *
 563  * If this is the first call for this frame, evaluates the child element
 564  * there should be exactly one in \p g.
 565  * The compiler has taken care that the child actually stores the evaluated
 566  * value in the value pointer of this element.
 567  * Assumes that \p g is persistent for the duration of the whole evaluation.
 568  *
 569  * This function is used as gmx::SelectionTreeElement::evaluate for
 570  * \ref SEL_SUBEXPR elements that have a static evaluation group, and hence do
 571  * not need full subexpression handling.
 572  */
 573 void
 574 _gmx_sel_evaluate_subexpr_staticeval(gmx_sel_evaluate_t                *data,
 575                                      const SelectionTreeElementPointer &sel,
 576                                      gmx_ana_index_t                   *g)
 577 {
 578     if (sel->u.cgrp.isize == 0)
 579     {
 580         sel->child->evaluate(data, sel->child, g);
 581         sel->v.nr = sel->child->v.nr;
 582         if (!g)
 583         {
 584             sel->u.cgrp.isize = -1;
 585         }
 586         else
 587         {
 588             gmx_ana_index_set(&sel->u.cgrp, g->isize, g->index, sel->u.cgrp.name, 0);
 589         }
 590     }
 591 }
 592
 593 /*!
 594  * \param[in]  data  Data for the current frame.
 595  * \param[in]  sel   Selection element being evaluated.
 596  * \param[in]  g     Group for which \p sel should be evaluated.
 597  * \returns    0 on success, a non-zero error code on error.
 598  *
 599  * Finds the part of \p g for which the subexpression
 600  * has not yet been evaluated by comparing \p g to \p sel->u.cgrp.
 601  * If the part is not empty, the child expression is evaluated for this
 602  * part, and the results merged to the old values of the child.
 603  * The value of \p sel itself is undefined after the call.
 604  *
 605  * \todo
 606  * The call to gmx_ana_index_difference() can take quite a lot of unnecessary
 607  * time if the subexpression is evaluated either several times for the same
 608  * group or for completely distinct groups.
 609  * However, in the majority of cases, these situations occur when
 610  * _gmx_sel_evaluate_subexpr_staticeval() can be used, so this should not be a
 611  * major problem.
 612  */
 613 void
 614 _gmx_sel_evaluate_subexpr(gmx_sel_evaluate_t                *data,
 615                           const SelectionTreeElementPointer &sel,
 616                           gmx_ana_index_t                   *g)
 617 {
 618     gmx_ana_index_t      gmiss;
 619
 620     MempoolGroupReserver gmissreserver(data->mp);
 621     if (sel->u.cgrp.isize == 0)
 622     {
 623         {
 624             SelelemTemporaryValueAssigner assigner(sel->child, *sel);
 625             sel->child->evaluate(data, sel->child, g);
 626         }
 627         /* We need to keep the name for the cgrp across the copy to avoid
 628          * problems if g has a name set. */
 629         char *name = sel->u.cgrp.name;
 630         gmx_ana_index_copy(&sel->u.cgrp, g, false);
 631         sel->u.cgrp.name = name;
 632         gmiss.isize      = 0;
 633     }
 634     else
 635     {
 636         gmissreserver.reserve(&gmiss, g->isize);
 637         gmx_ana_index_difference(&gmiss, g, &sel->u.cgrp);
 638         gmiss.name = NULL;
 639     }
 640     if (gmiss.isize > 0)
 641     {
 642         MempoolSelelemReserver reserver(sel->child, gmiss.isize);
 643         /* Evaluate the missing values for the child */
 644         sel->child->evaluate(data, sel->child, &gmiss);
 645         /* Merge the missing values to the existing ones. */
 646         if (sel->v.type == GROUP_VALUE)
 647         {
 648             gmx_ana_index_merge(sel->v.u.g, sel->child->v.u.g, sel->v.u.g);
 649         }
 650         else
 651         {
 652             int  i, j, k;
 653
 654             i = sel->u.cgrp.isize - 1;
 655             j = gmiss.isize - 1;
 656             /* TODO: This switch is kind of ugly, but it may be difficult to
 657              * do this portably without C++ templates. */
 658             switch (sel->v.type)
 659             {
 660                 case INT_VALUE:
 661                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 662                     {
 663                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 664                         {
 665                             sel->v.u.i[k] = sel->child->v.u.i[j--];
 666                         }
 667                         else
 668                         {
 669                             sel->v.u.i[k] = sel->v.u.i[i--];
 670                         }
 671                     }
 672                     break;
 673
 674                 case REAL_VALUE:
 675                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 676                     {
 677                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 678                         {
 679                             sel->v.u.r[k] = sel->child->v.u.r[j--];
 680                         }
 681                         else
 682                         {
 683                             sel->v.u.r[k] = sel->v.u.r[i--];
 684                         }
 685                     }
 686                     break;
 687
 688                 case STR_VALUE:
 689                     // Note: with the currently allowed syntax, this case is never
 690                     // reached.
 691                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 692                     {
 693                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 694                         {
 695                             sel->v.u.s[k] = sel->child->v.u.s[j--];
 696                         }
 697                         else
 698                         {
 699                             sel->v.u.s[k] = sel->v.u.s[i--];
 700                         }
 701                     }
 702                     break;
 703
 704                 case POS_VALUE:
 705                     /* TODO: Implement this */
 706                     GMX_THROW(gmx::NotImplementedError("position subexpressions not implemented properly"));
 707
 708                 case NO_VALUE:
 709                 case GROUP_VALUE:
 710                     GMX_THROW(gmx::InternalError("Invalid subexpression type"));
 711             }
 712         }
 713         gmx_ana_index_merge(&sel->u.cgrp, &sel->u.cgrp, &gmiss);
 714     }
 715 }
 716
 717 /*!
 718  * \param[in] data Data for the current frame.
 719  * \param[in] sel Selection element being evaluated.
 720  * \param[in] g   Group for which \p sel should be evaluated.
 721  * \returns   0 for success.
 722  *
 723  * Sets the value pointers of the child and its child to point to the same
 724  * memory as the value pointer of this element to avoid copying, and then
 725  * evaluates evaluates the child.
 726  *
 727  * This function is used as gmx::SelectionTreeElement:evaluate for
 728  * \ref SEL_SUBEXPRREF elements for which the \ref SEL_SUBEXPR does not have
 729  * other references.
 730  */
 731 void
 732 _gmx_sel_evaluate_subexprref_simple(gmx_sel_evaluate_t                *data,
 733                                     const SelectionTreeElementPointer &sel,
 734                                     gmx_ana_index_t                   *g)
 735 {
 736     if (g)
 737     {
 738         _gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
 739         _gmx_selvalue_setstore_alloc(&sel->child->child->v, sel->v.u.ptr,
 740                                      sel->child->child->v.nalloc);
 741         sel->child->evaluate(data, sel->child, g);
 742     }
 743     sel->v.nr = sel->child->v.nr;
 744     if (sel->u.param)
 745     {
 746         sel->u.param->val.nr = sel->v.nr;
 747         if (sel->u.param->nvalptr)
 748         {
 749             *sel->u.param->nvalptr = sel->u.param->val.nr;
 750         }
 751     }
 752 }
 753
 754 /*!
 755  * \param[in] data Data for the current frame.
 756  * \param[in] sel Selection element being evaluated.
 757  * \param[in] g   Group for which \p sel should be evaluated.
 758  * \returns   0 on success, a non-zero error code on error.
 759  *
 760  * If the value type is \ref POS_VALUE, the value of the child is simply
 761  * copied to set the value of \p sel (the child subexpression should
 762  * already have been evaluated by its root).
 763  * If the value type is something else, the child is evaluated for the
 764  * group \p g, and the value of the child is then copied.
 765  * There should be only one child element.
 766  *
 767  * This function is used as gmx::SelectionTreeElement::evaluate for
 768  * \ref SEL_SUBEXPRREF elements.
 769  */
 770 void
 771 _gmx_sel_evaluate_subexprref(gmx_sel_evaluate_t                *data,
 772                              const SelectionTreeElementPointer &sel,
 773                              gmx_ana_index_t                   *g)
 774 {
 775     int        i, j;
 776
 777     if (g)
 778     {
 779         sel->child->evaluate(data, sel->child, g);
 780     }
 781     const SelectionTreeElementPointer &expr = sel->child;
 782     switch (sel->v.type)
 783     {
 784         case INT_VALUE:
 785             if (!g)
 786             {
 787                 sel->v.nr = expr->v.nr;
 788                 memcpy(sel->v.u.i, expr->v.u.i, sel->v.nr*sizeof(*sel->v.u.i));
 789             }
 790             else
 791             {
 792                 sel->v.nr = g->isize;
 793                 /* Extract the values corresponding to g */
 794                 for (i = j = 0; i < g->isize; ++i, ++j)
 795                 {
 796                     while (sel->child->u.cgrp.index[j] < g->index[i])
 797                     {
 798                         ++j;
 799                     }
 800                     sel->v.u.i[i] = expr->v.u.i[j];
 801                 }
 802             }
 803             break;
 804
 805         case REAL_VALUE:
 806             if (!g)
 807             {
 808                 sel->v.nr = expr->v.nr;
 809                 memcpy(sel->v.u.r, expr->v.u.r, sel->v.nr*sizeof(*sel->v.u.r));
 810             }
 811             else
 812             {
 813                 sel->v.nr = g->isize;
 814                 /* Extract the values corresponding to g */
 815                 for (i = j = 0; i < g->isize; ++i, ++j)
 816                 {
 817                     while (sel->child->u.cgrp.index[j] < g->index[i])
 818                     {
 819                         ++j;
 820                     }
 821                     sel->v.u.r[i] = expr->v.u.r[j];
 822                 }
 823             }
 824             break;
 825
 826         case STR_VALUE:
 827             if (!g)
 828             {
 829                 sel->v.nr = expr->v.nr;
 830                 memcpy(sel->v.u.s, expr->v.u.s, sel->v.nr*sizeof(*sel->v.u.s));
 831             }
 832             else
 833             {
 834                 sel->v.nr = g->isize;
 835                 /* Extract the values corresponding to g */
 836                 for (i = j = 0; i < g->isize; ++i, ++j)
 837                 {
 838                     while (sel->child->u.cgrp.index[j] < g->index[i])
 839                     {
 840                         ++j;
 841                     }
 842                     sel->v.u.s[i] = expr->v.u.s[j];
 843                 }
 844             }
 845             break;
 846
 847         case POS_VALUE:
 848             /* Currently, there is no need to do anything fancy here,
 849              * but some future extensions may need a more flexible
 850              * implementation. */
 851             gmx_ana_pos_copy(sel->v.u.p, expr->v.u.p, false);
 852             break;
 853
 854         case GROUP_VALUE:
 855             if (!g)
 856             {
 857                 gmx_ana_index_copy(sel->v.u.g, expr->v.u.g, false);
 858             }
 859             else
 860             {
 861                 gmx_ana_index_intersection(sel->v.u.g, expr->v.u.g, g);
 862             }
 863             break;
 864
 865         default: /* should not be reached */
 866             GMX_THROW(gmx::InternalError("Invalid subexpression reference type"));
 867     }
 868     /* Store the number of values if needed */
 869     if (sel->u.param)
 870     {
 871         sel->u.param->val.nr = sel->v.nr;
 872         if (sel->u.param->nvalptr)
 873         {
 874             *sel->u.param->nvalptr = sel->u.param->val.nr;
 875         }
 876     }
 877 }
 878
 879 /********************************************************************
 880  * METHOD EXPRESSION EVALUATION
 881  ********************************************************************/
 882
 883 /*!
 884  * \param[in] data Data for the current frame.
 885  * \param[in] sel Selection element being evaluated.
 886  * \param[in] g   Group for which \p sel should be evaluated.
 887  * \returns   0 on success, a non-zero error code on error.
 888  *
 889  * Evaluates each child of a \ref SEL_EXPRESSION element.
 890  * The value of \p sel is not touched.
 891  *
 892  * This function is not used as gmx::SelectionTreeElement::evaluate,
 893  * but is used internally.
 894  */
 895 void
 896 _gmx_sel_evaluate_method_params(gmx_sel_evaluate_t                *data,
 897                                 const SelectionTreeElementPointer &sel,
 898                                 gmx_ana_index_t                   *g)
 899 {
 900     SelectionTreeElementPointer child = sel->child;
 901     while (child)
 902     {
 903         if (child->evaluate && !(child->flags & SEL_EVALFRAME))
 904         {
 905             if (child->flags & SEL_ATOMVAL)
 906             {
 907                 child->evaluate(data, child, g);
 908             }
 909             else
 910             {
 911                 child->flags |= SEL_EVALFRAME;
 912                 child->evaluate(data, child, NULL);
 913             }
 914         }
 915         child = child->next;
 916     }
 917 }
 918
 919 /*!
 920  * \param[in] data Data for the current frame.
 921  * \param[in] sel Selection element being evaluated.
 922  * \param[in] g   Group for which \p sel should be evaluated.
 923  * \returns   0 on success, a non-zero error code on error.
 924  *
 925  * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
 926  * to evaluate any parameter values.
 927  * If this is the first time this expression is evaluated for
 928  * the frame, sel_framefunc() callback is called if one is provided.
 929  * If a reference position calculation has been initialized for this element,
 930  * the positions are also updated, and sel_updatefunc_pos() is used to
 931  * evaluate the value. Otherwise, sel_updatefunc() is used.
 932  *
 933  * This function is used as gmx::SelectionTreeElement::evaluate for
 934  * \ref SEL_EXPRESSION elements.
 935  */
 936 void
 937 _gmx_sel_evaluate_method(gmx_sel_evaluate_t                *data,
 938                          const SelectionTreeElementPointer &sel,
 939                          gmx_ana_index_t                   *g)
 940 {
 941     _gmx_sel_evaluate_method_params(data, sel, g);
 942     if (sel->flags & SEL_INITFRAME)
 943     {
 944         sel->flags &= ~SEL_INITFRAME;
 945         sel->u.expr.method->init_frame(data->top, data->fr, data->pbc,
 946                                        sel->u.expr.mdata);
 947     }
 948     if (sel->u.expr.pc)
 949     {
 950         gmx_ana_poscalc_update(sel->u.expr.pc, sel->u.expr.pos, g,
 951                                data->fr, data->pbc);
 952         sel->u.expr.method->pupdate(data->top, data->fr, data->pbc,
 953                                     sel->u.expr.pos, &sel->v,
 954                                     sel->u.expr.mdata);
 955     }
 956     else
 957     {
 958         sel->u.expr.method->update(data->top, data->fr, data->pbc, g,
 959                                    &sel->v, sel->u.expr.mdata);
 960     }
 961 }
 962
 963 /*!
 964  * \param[in] data Data for the current frame.
 965  * \param[in] sel Selection element being evaluated.
 966  * \param[in] g   Group for which \p sel should be evaluated.
 967  * \returns   0 on success, a non-zero error code on error.
 968  *
 969  * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
 970  * to evaluate any parameter values.
 971  * If this is the first time this expression is evaluated for
 972  * the frame, sel_framefunc() callback is called if one is provided.
 973  * The modifier is then evaluated using sel_updatefunc_pos().
 974  *
 975  * This function is used as gmx::SelectionTreeElement::evaluate for
 976  * \ref SEL_MODIFIER elements.
 977  */
 978 void
 979 _gmx_sel_evaluate_modifier(gmx_sel_evaluate_t                *data,
 980                            const SelectionTreeElementPointer &sel,
 981                            gmx_ana_index_t                   *g)
 982 {
 983     _gmx_sel_evaluate_method_params(data, sel, g);
 984     if (sel->flags & SEL_INITFRAME)
 985     {
 986         sel->flags &= ~SEL_INITFRAME;
 987         sel->u.expr.method->init_frame(data->top, data->fr, data->pbc,
 988                                        sel->u.expr.mdata);
 989     }
 990     GMX_RELEASE_ASSERT(sel->child != NULL,
 991                        "Modifier element with a value must have a child");
 992     if (sel->child->v.type != POS_VALUE)
 993     {
 994         GMX_THROW(gmx::NotImplementedError("Non-position valued modifiers not implemented"));
 995     }
 996     sel->u.expr.method->pupdate(data->top, data->fr, data->pbc,
 997                                 sel->child->v.u.p,
 998                                 &sel->v, sel->u.expr.mdata);
 999 }
1000
1001
1002 /********************************************************************
1003  * BOOLEAN EXPRESSION EVALUATION
1004  ********************************************************************/
1005
1006 /*!
1007  * \param[in] data Data for the current frame.
1008  * \param[in] sel Selection element being evaluated.
1009  * \param[in] g   Group for which \p sel should be evaluated.
1010  * \returns   0 on success, a non-zero error code on error.
1011  *
1012  * Evaluates the child element (there should be only one) in the group
1013  * \p g, and then sets the value of \p sel to the complement of the
1014  * child value.
1015  *
1016  * This function is used as gmx::SelectionTreeElement::evaluate for
1017  * \ref SEL_BOOLEAN elements with \ref BOOL_NOT.
1018  */
1019 void
1020 _gmx_sel_evaluate_not(gmx_sel_evaluate_t                *data,
1021                       const SelectionTreeElementPointer &sel,
1022                       gmx_ana_index_t                   *g)
1023 {
1024     MempoolSelelemReserver reserver(sel->child, g->isize);
1025     sel->child->evaluate(data, sel->child, g);
1026     gmx_ana_index_difference(sel->v.u.g, g, sel->child->v.u.g);
1027 }
1028
1029 /*!
1030  * \param[in] data Data for the current frame.
1031  * \param[in] sel Selection element being evaluated.
1032  * \param[in] g   Group for which \p sel should be evaluated.
1033  * \returns   0 on success, a non-zero error code on error.
1034  *
1035  * Short-circuiting evaluation of logical AND expressions.
1036  *
1037  * Starts by evaluating the first child element in the group \p g.
1038  * The each following child element is evaluated in the intersection
1039  * of all the previous values until all children have been evaluated
1040  * or the intersection becomes empty.
1041  * The value of \p sel is set to the intersection of all the (evaluated)
1042  * child values.
1043  *
1044  * If the first child does not have an evaluation function, it is skipped
1045  * and the evaluation is started at the second child.
1046  * This happens if the first child is a constant expression and during
1047  * compilation it was detected that the evaluation group is always a subset
1048  * of the constant group
1049  * (currently, the compiler never detects this).
1050  *
1051  * This function is used as gmx::SelectionTreeElement::evaluate for
1052  * \ref SEL_BOOLEAN elements with \ref BOOL_AND.
1053  */
1054 void
1055 _gmx_sel_evaluate_and(gmx_sel_evaluate_t                *data,
1056                       const SelectionTreeElementPointer &sel,
1057                       gmx_ana_index_t                   *g)
1058 {
1059     SelectionTreeElementPointer child = sel->child;
1060     /* Skip the first child if it does not have an evaluation function. */
1061     if (!child->evaluate)
1062     {
1063         child = child->next;
1064     }
1065     /* Evaluate the first child */
1066     {
1067         MempoolSelelemReserver reserver(child, g->isize);
1068         child->evaluate(data, child, g);
1069         gmx_ana_index_copy(sel->v.u.g, child->v.u.g, false);
1070     }
1071     child = child->next;
1072     while (child && sel->v.u.g->isize > 0)
1073     {
1074         MempoolSelelemReserver reserver(child, sel->v.u.g->isize);
1075         child->evaluate(data, child, sel->v.u.g);
1076         gmx_ana_index_intersection(sel->v.u.g, sel->v.u.g, child->v.u.g);
1077         child = child->next;
1078     }
1079 }
1080
1081 /*!
1082  * \param[in] data Data for the current frame.
1083  * \param[in] sel Selection element being evaluated.
1084  * \param[in] g   Group for which \p sel should be evaluated.
1085  * \returns   0 on success, a non-zero error code on error.
1086  *
1087  * Short-circuiting evaluation of logical OR expressions.
1088  *
1089  * Starts by evaluating the first child element in the group \p g.
1090  * For each subsequent child, finds the part of \p g that is not
1091  * included the value of any previous child, and evaluates the child
1092  * in that group until the last child is evaluated or all of \p g
1093  * is included in some child value.
1094  * The value of \p sel is set to the union of all the (evaluated)
1095  * child values.
1096  *
1097  * If the first child does not have an evaluation function, its value is
1098  * used without evaluation.
1099  * This happens if the first child is a constant expression, the selection
1100  * has been compiled, and the evaluation group is the same for each frame.
1101  * In this case, the compiler has taken care of that the child value is a
1102  * subset of \p g, making it unnecessary to evaluate it.
1103  *
1104  * This function is used as gmx::SelectionTreeElement::evaluate for
1105  * \ref SEL_BOOLEAN elements with \ref BOOL_OR.
1106  */
1107 void
1108 _gmx_sel_evaluate_or(gmx_sel_evaluate_t                *data,
1109                      const SelectionTreeElementPointer &sel,
1110                      gmx_ana_index_t                   *g)
1111 {
1112     gmx_ana_index_t             tmp, tmp2;
1113
1114     SelectionTreeElementPointer child = sel->child;
1115     if (child->evaluate)
1116     {
1117         MempoolSelelemReserver reserver(child, g->isize);
1118         child->evaluate(data, child, g);
1119         gmx_ana_index_partition(sel->v.u.g, &tmp, g, child->v.u.g);
1120     }
1121     else
1122     {
1123         gmx_ana_index_partition(sel->v.u.g, &tmp, g, child->v.u.g);
1124     }
1125     child = child->next;
1126     while (child && tmp.isize > 0)
1127     {
1128         tmp.name = NULL;
1129         {
1130             MempoolSelelemReserver reserver(child, tmp.isize);
1131             child->evaluate(data, child, &tmp);
1132             gmx_ana_index_partition(&tmp, &tmp2, &tmp, child->v.u.g);
1133         }
1134         sel->v.u.g->isize += tmp.isize;
1135         tmp.isize          = tmp2.isize;
1136         tmp.index          = tmp2.index;
1137         child              = child->next;
1138     }
1139     gmx_ana_index_sort(sel->v.u.g);
1140 }
1141
1142
1143 /********************************************************************
1144  * ARITHMETIC EVALUATION
1145  ********************************************************************/
1146
1147 /*!
1148  * \param[in] data Data for the current frame.
1149  * \param[in] sel  Selection element being evaluated.
1150  * \param[in] g    Group for which \p sel should be evaluated.
1151  * \returns   0 on success, a non-zero error code on error.
1152  */
1153 void
1154 _gmx_sel_evaluate_arithmetic(gmx_sel_evaluate_t                *data,
1155                              const SelectionTreeElementPointer &sel,
1156                              gmx_ana_index_t                   *g)
1157 {
1158     int         n, i, i1, i2;
1159     real        lval, rval = 0., val = 0.;
1160
1161     const SelectionTreeElementPointer &left  = sel->child;
1162     const SelectionTreeElementPointer &right = left->next;
1163
1164     SelelemTemporaryValueAssigner      assigner;
1165     MempoolSelelemReserver             reserver;
1166     if (left->mempool)
1167     {
1168         assigner.assign(left, *sel);
1169         if (right)
1170         {
1171             reserver.reserve(right, g->isize);
1172         }
1173     }
1174     else if (right && right->mempool)
1175     {
1176         assigner.assign(right, *sel);
1177     }
1178     _gmx_sel_evaluate_children(data, sel, g);
1179
1180     n         = (sel->flags & SEL_SINGLEVAL) ? 1 : g->isize;
1181     sel->v.nr = n;
1182
1183     bool bArithNeg = (sel->u.arith.type == ARITH_NEG);
1184     GMX_ASSERT(right || bArithNeg,
1185                "Right operand cannot be null except for negations");
1186     for (i = i1 = i2 = 0; i < n; ++i)
1187     {
1188         lval = left->v.u.r[i1];
1189         if (!bArithNeg)
1190         {
1191             rval = right->v.u.r[i2];
1192         }
1193         switch (sel->u.arith.type)
1194         {
1195             case ARITH_PLUS:    val = lval + rval;     break;
1196             case ARITH_MINUS:   val = lval - rval;     break;
1197             case ARITH_NEG:     val = -lval;           break;
1198             case ARITH_MULT:    val = lval * rval;     break;
1199             case ARITH_DIV:     val = lval / rval;     break;
1200             case ARITH_EXP:     val = pow(lval, rval); break;
1201         }
1202         sel->v.u.r[i] = val;
1203         if (!(left->flags & SEL_SINGLEVAL))
1204         {
1205             ++i1;
1206         }
1207         if (!bArithNeg && !(right->flags & SEL_SINGLEVAL))
1208         {
1209             ++i2;
1210         }
1211     }
1212 }