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