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