src/gromacs/selection/evaluate.cpp

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