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