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