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 void SelectionEvaluator::evaluate(SelectionCollection* coll, t_trxframe* fr, t_pbc* pbc)
 423 {
 424     gmx_ana_selcollection_t* sc = &coll->impl_->sc_;
 425     gmx_sel_evaluate_t       data;
 426
 427     _gmx_sel_evaluate_init(&data, sc->mempool, &sc->gall, sc->top, fr, pbc);
 428     init_frame_eval(sc->root);
 429     SelectionTreeElementPointer sel = sc->root;
 430     while (sel)
 431     {
 432         /* Clear the evaluation group of subexpressions */
 433         if (sel->child && sel->child->type == SEL_SUBEXPR && sel->child->evaluate != nullptr)
 434         {
 435             sel->child->u.cgrp.isize = 0;
 436             /* Not strictly necessary, because the value will be overwritten
 437              * during first evaluation of the subexpression anyways, but we
 438              * clear the group for clarity. Note that this is _not_ done during
 439              * compilation because of some additional complexities involved
 440              * (see compiler.cpp), so it should not be relied upon in
 441              * _gmx_sel_evaluate_subexpr(). */
 442             if (sel->child->v.type == GROUP_VALUE)
 443             {
 444                 sel->child->v.u.g->isize = 0;
 445             }
 446         }
 447         if (sel->evaluate)
 448         {
 449             sel->evaluate(&data, sel, nullptr);
 450         }
 451         sel = sel->next;
 452     }
 453     /* Update selection information */
 454     SelectionDataList::const_iterator isel;
 455     for (isel = sc->sel.begin(); isel != sc->sel.end(); ++isel)
 456     {
 457         internal::SelectionData& sel = **isel;
 458         sel.refreshMassesAndCharges(sc->top);
 459         sel.updateCoveredFractionForFrame();
 460     }
 461 }
 462
 463 /*!
 464  * \param[in,out] coll  The selection collection to evaluate.
 465  * \param[in]     nframes Total number of frames.
 466  */
 467 void SelectionEvaluator::evaluateFinal(SelectionCollection* coll, int nframes)
 468 {
 469     gmx_ana_selcollection_t* sc = &coll->impl_->sc_;
 470
 471     SelectionDataList::const_iterator isel;
 472     for (isel = sc->sel.begin(); isel != sc->sel.end(); ++isel)
 473     {
 474         internal::SelectionData& sel = **isel;
 475         sel.restoreOriginalPositions(sc->top);
 476         sel.computeAverageCoveredFraction(nframes);
 477     }
 478 }
 479
 480 } // namespace gmx
 481
 482 /*!
 483  * \param[in] data Data for the current frame.
 484  * \param[in] sel  Selection element being evaluated.
 485  * \param[in] g    Group for which \p sel should be evaluated.
 486  * \returns   0 on success, a non-zero error code on error.
 487  *
 488  * Evaluates each child of \p sel in \p g.
 489  */
 490 void _gmx_sel_evaluate_children(gmx_sel_evaluate_t*                     data,
 491                                 const gmx::SelectionTreeElementPointer& sel,
 492                                 gmx_ana_index_t*                        g)
 493 {
 494     SelectionTreeElementPointer child = sel->child;
 495     while (child)
 496     {
 497         if (child->evaluate)
 498         {
 499             child->evaluate(data, child, g);
 500         }
 501         child = child->next;
 502     }
 503 }
 504
 505 void _gmx_sel_evaluate_root(gmx_sel_evaluate_t*                     data,
 506                             const gmx::SelectionTreeElementPointer& sel,
 507                             gmx_ana_index_t* /* g */)
 508 {
 509     if (sel->u.cgrp.isize == 0 || !sel->child->evaluate)
 510     {
 511         return;
 512     }
 513
 514     sel->child->evaluate(data, sel->child, sel->u.cgrp.isize < 0 ? nullptr : &sel->u.cgrp);
 515 }
 516
 517 void _gmx_sel_evaluate_static(gmx_sel_evaluate_t* /* data */,
 518                               const gmx::SelectionTreeElementPointer& sel,
 519                               gmx_ana_index_t*                        g)
 520 {
 521     if (sel->flags & SEL_UNSORTED)
 522     {
 523         // This only works if g contains all the atoms, but that is currently
 524         // the only supported case.
 525         gmx_ana_index_copy(sel->v.u.g, &sel->u.cgrp, false);
 526     }
 527     else
 528     {
 529         gmx_ana_index_intersection(sel->v.u.g, &sel->u.cgrp, g);
 530     }
 531 }
 532
 533
 534 /*********************************************************************
 535  * SUBEXPRESSION EVALUATION
 536  *********************************************************************/
 537
 538 /*!
 539  * \param[in] data Data for the current frame.
 540  * \param[in] sel  Selection element being evaluated.
 541  * \param[in] g    Group for which \p sel should be evaluated.
 542  * \returns   0 on success, a non-zero error code on error.
 543  *
 544  * Evaluates the child element (there should be exactly one) in \p g.
 545  * The compiler has taken care that the child actually stores the evaluated
 546  * value in the value pointer of this element.
 547  *
 548  * This function is used as gmx::SelectionTreeElement::evaluate for
 549  * \ref SEL_SUBEXPR elements that are used only once, and hence do not need
 550  * full subexpression handling.
 551  */
 552 void _gmx_sel_evaluate_subexpr_simple(gmx_sel_evaluate_t*                     data,
 553                                       const gmx::SelectionTreeElementPointer& sel,
 554                                       gmx_ana_index_t*                        g)
 555 {
 556     if (sel->child->evaluate)
 557     {
 558         sel->child->evaluate(data, sel->child, g);
 559     }
 560     sel->v.nr = sel->child->v.nr;
 561 }
 562
 563 /*!
 564  * \param[in] data Data for the current frame.
 565  * \param[in] sel  Selection element being evaluated.
 566  * \param[in] g    Group for which \p sel should be evaluated.
 567  * \returns   0 on success, a non-zero error code on error.
 568  *
 569  * If this is the first call for this frame, evaluates the child element
 570  * there should be exactly one in \p g.
 571  * The compiler has taken care that the child actually stores the evaluated
 572  * value in the value pointer of this element.
 573  * Assumes that \p g is persistent for the duration of the whole evaluation.
 574  *
 575  * This function is used as gmx::SelectionTreeElement::evaluate for
 576  * \ref SEL_SUBEXPR elements that have a static evaluation group, and hence do
 577  * not need full subexpression handling.
 578  */
 579 void _gmx_sel_evaluate_subexpr_staticeval(gmx_sel_evaluate_t*                     data,
 580                                           const gmx::SelectionTreeElementPointer& sel,
 581                                           gmx_ana_index_t*                        g)
 582 {
 583     if (sel->u.cgrp.isize == 0)
 584     {
 585         sel->child->evaluate(data, sel->child, g);
 586         sel->v.nr = sel->child->v.nr;
 587         if (!g)
 588         {
 589             sel->u.cgrp.isize = -1;
 590         }
 591         else
 592         {
 593             gmx_ana_index_set(&sel->u.cgrp, g->isize, g->index, 0);
 594         }
 595     }
 596 }
 597
 598 /*!
 599  * \param[in]  data  Data for the current frame.
 600  * \param[in]  sel   Selection element being evaluated.
 601  * \param[in]  g     Group for which \p sel should be evaluated.
 602  * \returns    0 on success, a non-zero error code on error.
 603  *
 604  * Finds the part of \p g for which the subexpression
 605  * has not yet been evaluated by comparing \p g to \p sel->u.cgrp.
 606  * If the part is not empty, the child expression is evaluated for this
 607  * part, and the results merged to the old values of the child.
 608  * The value of \p sel itself is undefined after the call.
 609  *
 610  * \todo
 611  * The call to gmx_ana_index_difference() can take quite a lot of unnecessary
 612  * time if the subexpression is evaluated either several times for the same
 613  * group or for completely distinct groups.
 614  * However, in the majority of cases, these situations occur when
 615  * _gmx_sel_evaluate_subexpr_staticeval() can be used, so this should not be a
 616  * major problem.
 617  */
 618 void _gmx_sel_evaluate_subexpr(gmx_sel_evaluate_t*                     data,
 619                                const gmx::SelectionTreeElementPointer& sel,
 620                                gmx_ana_index_t*                        g)
 621 {
 622     gmx_ana_index_t gmiss;
 623
 624     MempoolGroupReserver gmissreserver(data->mp);
 625     if (sel->u.cgrp.isize == 0)
 626     {
 627         {
 628             SelelemTemporaryValueAssigner assigner(sel->child, *sel);
 629             sel->child->evaluate(data, sel->child, g);
 630         }
 631         gmx_ana_index_copy(&sel->u.cgrp, g, false);
 632         gmiss.isize = 0;
 633     }
 634     else
 635     {
 636         gmissreserver.reserve(&gmiss, g->isize);
 637         gmx_ana_index_difference(&gmiss, g, &sel->u.cgrp);
 638     }
 639     if (gmiss.isize > 0)
 640     {
 641         MempoolSelelemReserver reserver(sel->child, gmiss.isize);
 642         /* Evaluate the missing values for the child */
 643         sel->child->evaluate(data, sel->child, &gmiss);
 644         /* Merge the missing values to the existing ones. */
 645         if (sel->v.type == GROUP_VALUE)
 646         {
 647             gmx_ana_index_merge(sel->v.u.g, sel->child->v.u.g, sel->v.u.g);
 648         }
 649         else
 650         {
 651             int i, j, k;
 652
 653             i = sel->u.cgrp.isize - 1;
 654             j = gmiss.isize - 1;
 655             /* TODO: This switch is kind of ugly, but it may be difficult to
 656              * do this portably without C++ templates. */
 657             switch (sel->v.type)
 658             {
 659                 case INT_VALUE:
 660                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 661                     {
 662                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 663                         {
 664                             sel->v.u.i[k] = sel->child->v.u.i[j--];
 665                         }
 666                         else
 667                         {
 668                             sel->v.u.i[k] = sel->v.u.i[i--];
 669                         }
 670                     }
 671                     break;
 672
 673                 case REAL_VALUE:
 674                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 675                     {
 676                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 677                         {
 678                             sel->v.u.r[k] = sel->child->v.u.r[j--];
 679                         }
 680                         else
 681                         {
 682                             sel->v.u.r[k] = sel->v.u.r[i--];
 683                         }
 684                     }
 685                     break;
 686
 687                 case STR_VALUE:
 688                     // Note: with the currently allowed syntax, this case is never
 689                     // reached.
 690                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 691                     {
 692                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 693                         {
 694                             sel->v.u.s[k] = sel->child->v.u.s[j--];
 695                         }
 696                         else
 697                         {
 698                             sel->v.u.s[k] = sel->v.u.s[i--];
 699                         }
 700                     }
 701                     break;
 702
 703                 case POS_VALUE:
 704                     /* TODO: Implement this */
 705                     GMX_THROW(gmx::NotImplementedError(
 706                             "position subexpressions not implemented properly"));
 707
 708                 case NO_VALUE:
 709                 case GROUP_VALUE: GMX_THROW(gmx::InternalError("Invalid subexpression type"));
 710             }
 711         }
 712         gmx_ana_index_merge(&sel->u.cgrp, &sel->u.cgrp, &gmiss);
 713     }
 714 }
 715
 716 /*!
 717  * \param[in] data Data for the current frame.
 718  * \param[in] sel Selection element being evaluated.
 719  * \param[in] g   Group for which \p sel should be evaluated.
 720  * \returns   0 for success.
 721  *
 722  * Sets the value pointers of the child and its child to point to the same
 723  * memory as the value pointer of this element to avoid copying, and then
 724  * evaluates evaluates the child.
 725  *
 726  * This function is used as gmx::SelectionTreeElement:evaluate for
 727  * \ref SEL_SUBEXPRREF elements for which the \ref SEL_SUBEXPR does not have
 728  * other references.
 729  */
 730 void _gmx_sel_evaluate_subexprref_simple(gmx_sel_evaluate_t*                     data,
 731                                          const gmx::SelectionTreeElementPointer& sel,
 732                                          gmx_ana_index_t*                        g)
 733 {
 734     if (g)
 735     {
 736         _gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
 737         _gmx_selvalue_setstore_alloc(&sel->child->child->v, sel->v.u.ptr, sel->child->child->v.nalloc);
 738         sel->child->evaluate(data, sel->child, g);
 739     }
 740     sel->v.nr = sel->child->v.nr;
 741     if (sel->u.param)
 742     {
 743         sel->u.param->val.nr = sel->v.nr;
 744         if (sel->u.param->nvalptr)
 745         {
 746             *sel->u.param->nvalptr = sel->u.param->val.nr;
 747         }
 748     }
 749 }
 750
 751 /*!
 752  * \param[in] data Data for the current frame.
 753  * \param[in] sel Selection element being evaluated.
 754  * \param[in] g   Group for which \p sel should be evaluated.
 755  * \returns   0 on success, a non-zero error code on error.
 756  *
 757  * If the value type is \ref POS_VALUE, the value of the child is simply
 758  * copied to set the value of \p sel (the child subexpression should
 759  * already have been evaluated by its root).
 760  * If the value type is something else, the child is evaluated for the
 761  * group \p g, and the value of the child is then copied.
 762  * There should be only one child element.
 763  *
 764  * This function is used as gmx::SelectionTreeElement::evaluate for
 765  * \ref SEL_SUBEXPRREF elements.
 766  */
 767 void _gmx_sel_evaluate_subexprref(gmx_sel_evaluate_t*                     data,
 768                                   const gmx::SelectionTreeElementPointer& sel,
 769                                   gmx_ana_index_t*                        g)
 770 {
 771     int i, j;
 772
 773     if (g != nullptr && sel->child->evaluate != nullptr)
 774     {
 775         sel->child->evaluate(data, sel->child, g);
 776     }
 777     const SelectionTreeElementPointer& expr = sel->child;
 778     switch (sel->v.type)
 779     {
 780         case INT_VALUE:
 781             if (!g)
 782             {
 783                 sel->v.nr = expr->v.nr;
 784                 memcpy(sel->v.u.i, expr->v.u.i, sel->v.nr * sizeof(*sel->v.u.i));
 785             }
 786             else
 787             {
 788                 sel->v.nr = g->isize;
 789                 /* Extract the values corresponding to g */
 790                 for (i = j = 0; i < g->isize; ++i, ++j)
 791                 {
 792                     while (sel->child->u.cgrp.index[j] < g->index[i])
 793                     {
 794                         ++j;
 795                     }
 796                     sel->v.u.i[i] = expr->v.u.i[j];
 797                 }
 798             }
 799             break;
 800
 801         case REAL_VALUE:
 802             if (!g)
 803             {
 804                 sel->v.nr = expr->v.nr;
 805                 memcpy(sel->v.u.r, expr->v.u.r, sel->v.nr * sizeof(*sel->v.u.r));
 806             }
 807             else
 808             {
 809                 sel->v.nr = g->isize;
 810                 /* Extract the values corresponding to g */
 811                 for (i = j = 0; i < g->isize; ++i, ++j)
 812                 {
 813                     while (sel->child->u.cgrp.index[j] < g->index[i])
 814                     {
 815                         ++j;
 816                     }
 817                     sel->v.u.r[i] = expr->v.u.r[j];
 818                 }
 819             }
 820             break;
 821
 822         case STR_VALUE:
 823             if (!g)
 824             {
 825                 sel->v.nr = expr->v.nr;
 826                 memcpy(sel->v.u.s, expr->v.u.s, sel->v.nr * sizeof(*sel->v.u.s));
 827             }
 828             else
 829             {
 830                 sel->v.nr = g->isize;
 831                 /* Extract the values corresponding to g */
 832                 for (i = j = 0; i < g->isize; ++i, ++j)
 833                 {
 834                     while (sel->child->u.cgrp.index[j] < g->index[i])
 835                     {
 836                         ++j;
 837                     }
 838                     sel->v.u.s[i] = expr->v.u.s[j];
 839                 }
 840             }
 841             break;
 842
 843         case POS_VALUE:
 844             /* Currently, there is no need to do anything fancy here,
 845              * but some future extensions may need a more flexible
 846              * implementation. */
 847             gmx_ana_pos_copy(sel->v.u.p, expr->v.u.p, false);
 848             break;
 849
 850         case GROUP_VALUE:
 851             if (!g)
 852             {
 853                 gmx_ana_index_copy(sel->v.u.g, expr->v.u.g, false);
 854             }
 855             else
 856             {
 857                 gmx_ana_index_intersection(sel->v.u.g, expr->v.u.g, g);
 858             }
 859             break;
 860
 861         default: /* should not be reached */
 862             GMX_THROW(gmx::InternalError("Invalid subexpression reference type"));
 863     }
 864     /* Store the number of values if needed */
 865     if (sel->u.param)
 866     {
 867         sel->u.param->val.nr = sel->v.nr;
 868         if (sel->u.param->nvalptr)
 869         {
 870             *sel->u.param->nvalptr = sel->u.param->val.nr;
 871         }
 872     }
 873 }
 874
 875 /********************************************************************
 876  * METHOD EXPRESSION EVALUATION
 877  ********************************************************************/
 878
 879 /*!
 880  * \param[in] data Data for the current frame.
 881  * \param[in] sel Selection element being evaluated.
 882  * \param[in] g   Group for which \p sel should be evaluated.
 883  * \returns   0 on success, a non-zero error code on error.
 884  *
 885  * Evaluates each child of a \ref SEL_EXPRESSION element.
 886  * The value of \p sel is not touched.
 887  *
 888  * This function is not used as gmx::SelectionTreeElement::evaluate,
 889  * but is used internally.
 890  */
 891 void _gmx_sel_evaluate_method_params(gmx_sel_evaluate_t*                     data,
 892                                      const gmx::SelectionTreeElementPointer& sel,
 893                                      gmx_ana_index_t*                        g)
 894 {
 895     SelectionTreeElementPointer child = sel->child;
 896     while (child)
 897     {
 898         if (child->evaluate && !(child->flags & SEL_EVALFRAME))
 899         {
 900             if (child->flags & SEL_ATOMVAL)
 901             {
 902                 child->evaluate(data, child, g);
 903             }
 904             else
 905             {
 906                 child->flags |= SEL_EVALFRAME;
 907                 child->evaluate(data, child, nullptr);
 908             }
 909         }
 910         child = child->next;
 911     }
 912 }
 913
 914 /*!
 915  * \param[in] data Data for the current frame.
 916  * \param[in] sel Selection element being evaluated.
 917  * \param[in] g   Group for which \p sel should be evaluated.
 918  * \returns   0 on success, a non-zero error code on error.
 919  *
 920  * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
 921  * to evaluate any parameter values.
 922  * If this is the first time this expression is evaluated for
 923  * the frame, sel_framefunc() callback is called if one is provided.
 924  * If a reference position calculation has been initialized for this element,
 925  * the positions are also updated, and sel_updatefunc_pos() is used to
 926  * evaluate the value. Otherwise, sel_updatefunc() is used.
 927  *
 928  * This function is used as gmx::SelectionTreeElement::evaluate for
 929  * \ref SEL_EXPRESSION elements.
 930  */
 931 void _gmx_sel_evaluate_method(gmx_sel_evaluate_t*                     data,
 932                               const gmx::SelectionTreeElementPointer& sel,
 933                               gmx_ana_index_t*                        g)
 934 {
 935     _gmx_sel_evaluate_method_params(data, sel, g);
 936     gmx::SelMethodEvalContext context(data->top, data->fr, data->pbc);
 937     if (sel->flags & SEL_INITFRAME)
 938     {
 939         sel->flags &= ~SEL_INITFRAME;
 940         sel->u.expr.method->init_frame(context, sel->u.expr.mdata);
 941     }
 942     if (sel->u.expr.pc)
 943     {
 944         gmx_ana_poscalc_update(sel->u.expr.pc, sel->u.expr.pos, g, data->fr, data->pbc);
 945         sel->u.expr.method->pupdate(context, sel->u.expr.pos, &sel->v, sel->u.expr.mdata);
 946         if ((sel->flags & SEL_ATOMVAL) && sel->v.nr < g->isize)
 947         {
 948             switch (sel->v.type)
 949             {
 950                 case REAL_VALUE:
 951                     expandValueForPositions(sel->v.u.r, &sel->v.nr, sel->u.expr.pos);
 952                     break;
 953                 default:
 954                     GMX_RELEASE_ASSERT(false,
 955                                        "Unimplemented value type for position update method");
 956             }
 957         }
 958     }
 959     else
 960     {
 961         sel->u.expr.method->update(context, g, &sel->v, sel->u.expr.mdata);
 962     }
 963 }
 964
 965 /*!
 966  * \param[in] data Data for the current frame.
 967  * \param[in] sel Selection element being evaluated.
 968  * \param[in] g   Group for which \p sel should be evaluated.
 969  * \returns   0 on success, a non-zero error code on error.
 970  *
 971  * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
 972  * to evaluate any parameter values.
 973  * If this is the first time this expression is evaluated for
 974  * the frame, sel_framefunc() callback is called if one is provided.
 975  * The modifier is then evaluated using sel_updatefunc_pos().
 976  *
 977  * This function is used as gmx::SelectionTreeElement::evaluate for
 978  * \ref SEL_MODIFIER elements.
 979  */
 980 void _gmx_sel_evaluate_modifier(gmx_sel_evaluate_t*                     data,
 981                                 const gmx::SelectionTreeElementPointer& sel,
 982                                 gmx_ana_index_t*                        g)
 983 {
 984     _gmx_sel_evaluate_method_params(data, sel, g);
 985     gmx::SelMethodEvalContext context(data->top, data->fr, data->pbc);
 986     if (sel->flags & SEL_INITFRAME)
 987     {
 988         sel->flags &= ~SEL_INITFRAME;
 989         sel->u.expr.method->init_frame(context, sel->u.expr.mdata);
 990     }
 991     if (sel->child && sel->child->v.type != POS_VALUE)
 992     {
 993         GMX_THROW(gmx::NotImplementedError("Non-position valued modifiers not implemented"));
 994     }
 995     sel->u.expr.method->pupdate(context, nullptr, &sel->v, sel->u.expr.mdata);
 996 }
 997
 998
 999 /********************************************************************
1000  * BOOLEAN EXPRESSION EVALUATION
1001  ********************************************************************/
1002
1003 /*!
1004  * \param[in] data Data for the current frame.
1005  * \param[in] sel Selection element being evaluated.
1006  * \param[in] g   Group for which \p sel should be evaluated.
1007  * \returns   0 on success, a non-zero error code on error.
1008  *
1009  * Evaluates the child element (there should be only one) in the group
1010  * \p g, and then sets the value of \p sel to the complement of the
1011  * child value.
1012  *
1013  * This function is used as gmx::SelectionTreeElement::evaluate for
1014  * \ref SEL_BOOLEAN elements with \ref BOOL_NOT.
1015  */
1016 void _gmx_sel_evaluate_not(gmx_sel_evaluate_t*                     data,
1017                            const gmx::SelectionTreeElementPointer& sel,
1018                            gmx_ana_index_t*                        g)
1019 {
1020     MempoolSelelemReserver reserver(sel->child, g->isize);
1021     sel->child->evaluate(data, sel->child, g);
1022     gmx_ana_index_difference(sel->v.u.g, g, sel->child->v.u.g);
1023 }
1024
1025 /*!
1026  * \param[in] data Data for the current frame.
1027  * \param[in] sel Selection element being evaluated.
1028  * \param[in] g   Group for which \p sel should be evaluated.
1029  * \returns   0 on success, a non-zero error code on error.
1030  *
1031  * Short-circuiting evaluation of logical AND expressions.
1032  *
1033  * Starts by evaluating the first child element in the group \p g.
1034  * The each following child element is evaluated in the intersection
1035  * of all the previous values until all children have been evaluated
1036  * or the intersection becomes empty.
1037  * The value of \p sel is set to the intersection of all the (evaluated)
1038  * child values.
1039  *
1040  * If the first child does not have an evaluation function, it is skipped
1041  * and the evaluation is started at the second child.
1042  * This happens if the first child is a constant expression and during
1043  * compilation it was detected that the evaluation group is always a subset
1044  * of the constant group
1045  * (currently, the compiler never detects this).
1046  *
1047  * This function is used as gmx::SelectionTreeElement::evaluate for
1048  * \ref SEL_BOOLEAN elements with \ref BOOL_AND.
1049  */
1050 void _gmx_sel_evaluate_and(gmx_sel_evaluate_t*                     data,
1051                            const gmx::SelectionTreeElementPointer& sel,
1052                            gmx_ana_index_t*                        g)
1053 {
1054     SelectionTreeElementPointer child = sel->child;
1055     /* Skip the first child if it does not have an evaluation function. */
1056     if (!child->evaluate)
1057     {
1058         child = child->next;
1059     }
1060     /* Evaluate the first child */
1061     {
1062         MempoolSelelemReserver reserver(child, g->isize);
1063         child->evaluate(data, child, g);
1064         gmx_ana_index_copy(sel->v.u.g, child->v.u.g, false);
1065     }
1066     child = child->next;
1067     while (child && sel->v.u.g->isize > 0)
1068     {
1069         MempoolSelelemReserver reserver(child, sel->v.u.g->isize);
1070         child->evaluate(data, child, sel->v.u.g);
1071         gmx_ana_index_intersection(sel->v.u.g, sel->v.u.g, child->v.u.g);
1072         child = child->next;
1073     }
1074 }
1075
1076 /*!
1077  * \param[in] data Data for the current frame.
1078  * \param[in] sel Selection element being evaluated.
1079  * \param[in] g   Group for which \p sel should be evaluated.
1080  * \returns   0 on success, a non-zero error code on error.
1081  *
1082  * Short-circuiting evaluation of logical OR expressions.
1083  *
1084  * Starts by evaluating the first child element in the group \p g.
1085  * For each subsequent child, finds the part of \p g that is not
1086  * included the value of any previous child, and evaluates the child
1087  * in that group until the last child is evaluated or all of \p g
1088  * is included in some child value.
1089  * The value of \p sel is set to the union of all the (evaluated)
1090  * child values.
1091  *
1092  * If the first child does not have an evaluation function, its value is
1093  * used without evaluation.
1094  * This happens if the first child is a constant expression, the selection
1095  * has been compiled, and the evaluation group is the same for each frame.
1096  * In this case, the compiler has taken care of that the child value is a
1097  * subset of \p g, making it unnecessary to evaluate it.
1098  *
1099  * This function is used as gmx::SelectionTreeElement::evaluate for
1100  * \ref SEL_BOOLEAN elements with \ref BOOL_OR.
1101  */
1102 void _gmx_sel_evaluate_or(gmx_sel_evaluate_t*                     data,
1103                           const gmx::SelectionTreeElementPointer& sel,
1104                           gmx_ana_index_t*                        g)
1105 {
1106     gmx_ana_index_t tmp, tmp2;
1107
1108     SelectionTreeElementPointer child = sel->child;
1109     if (child->evaluate)
1110     {
1111         MempoolSelelemReserver reserver(child, g->isize);
1112         child->evaluate(data, child, g);
1113         gmx_ana_index_partition(sel->v.u.g, &tmp, g, child->v.u.g);
1114     }
1115     else
1116     {
1117         gmx_ana_index_partition(sel->v.u.g, &tmp, g, child->v.u.g);
1118     }
1119     child = child->next;
1120     while (child && tmp.isize > 0)
1121     {
1122         {
1123             MempoolSelelemReserver reserver(child, tmp.isize);
1124             child->evaluate(data, child, &tmp);
1125             gmx_ana_index_partition(&tmp, &tmp2, &tmp, child->v.u.g);
1126         }
1127         sel->v.u.g->isize += tmp.isize;
1128         tmp.isize = tmp2.isize;
1129         tmp.index = tmp2.index;
1130         child     = child->next;
1131     }
1132     gmx_ana_index_sort(sel->v.u.g);
1133 }
1134
1135
1136 /********************************************************************
1137  * ARITHMETIC EVALUATION
1138  ********************************************************************/
1139
1140 /*!
1141  * \param[in] data Data for the current frame.
1142  * \param[in] sel  Selection element being evaluated.
1143  * \param[in] g    Group for which \p sel should be evaluated.
1144  * \returns   0 on success, a non-zero error code on error.
1145  */
1146 void _gmx_sel_evaluate_arithmetic(gmx_sel_evaluate_t*                     data,
1147                                   const gmx::SelectionTreeElementPointer& sel,
1148                                   gmx_ana_index_t*                        g)
1149 {
1150     int  n, i, i1, i2;
1151     real lval, rval = 0., val = 0.;
1152
1153     const SelectionTreeElementPointer& left  = sel->child;
1154     const SelectionTreeElementPointer& right = left->next;
1155
1156     SelelemTemporaryValueAssigner assigner;
1157     MempoolSelelemReserver        reserver;
1158     if (left->mempool)
1159     {
1160         assigner.assign(left, *sel);
1161         if (right)
1162         {
1163             reserver.reserve(right, g->isize);
1164         }
1165     }
1166     else if (right && right->mempool)
1167     {
1168         assigner.assign(right, *sel);
1169     }
1170     _gmx_sel_evaluate_children(data, sel, g);
1171
1172     n         = (sel->flags & SEL_SINGLEVAL) ? 1 : g->isize;
1173     sel->v.nr = n;
1174
1175     bool bArithNeg = (sel->u.arith.type == ARITH_NEG);
1176     GMX_ASSERT(right || bArithNeg, "Right operand cannot be null except for negations");
1177     for (i = i1 = i2 = 0; i < n; ++i)
1178     {
1179         lval = left->v.u.r[i1];
1180         if (!bArithNeg)
1181         {
1182             rval = right->v.u.r[i2];
1183         }
1184         switch (sel->u.arith.type)
1185         {
1186             case ARITH_PLUS: val = lval + rval; break;
1187             case ARITH_MINUS: val = lval - rval; break;
1188             case ARITH_NEG: val = -lval; break;
1189             case ARITH_MULT: val = lval * rval; break;
1190             case ARITH_DIV: val = lval / rval; break;
1191             case ARITH_EXP: val = pow(lval, rval); break;
1192         }
1193         sel->v.u.r[i] = val;
1194         if (!(left->flags & SEL_SINGLEVAL))
1195         {
1196             ++i1;
1197         }
1198         if (!bArithNeg && !(right->flags & SEL_SINGLEVAL))
1199         {
1200             ++i2;
1201         }
1202     }
1203 }