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,2021, 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         gmx_ana_index_reserve(sel->v.u.g, sel->u.cgrp.isize);
 526         // This only works if g contains all the atoms, but that is currently
 527         // the only supported case.
 528         gmx_ana_index_copy(sel->v.u.g, &sel->u.cgrp, false);
 529     }
 530     else
 531     {
 532         gmx_ana_index_intersection(sel->v.u.g, &sel->u.cgrp, g);
 533     }
 534 }
 535
 536
 537 /*********************************************************************
 538  * SUBEXPRESSION EVALUATION
 539  *********************************************************************/
 540
 541 /*!
 542  * \param[in] data Data for the current frame.
 543  * \param[in] sel  Selection element being evaluated.
 544  * \param[in] g    Group for which \p sel should be evaluated.
 545  * \returns   0 on success, a non-zero error code on error.
 546  *
 547  * Evaluates the child element (there should be exactly one) in \p g.
 548  * The compiler has taken care that the child actually stores the evaluated
 549  * value in the value pointer of this element.
 550  *
 551  * This function is used as gmx::SelectionTreeElement::evaluate for
 552  * \ref SEL_SUBEXPR elements that are used only once, and hence do not need
 553  * full subexpression handling.
 554  */
 555 void _gmx_sel_evaluate_subexpr_simple(gmx_sel_evaluate_t*                     data,
 556                                       const gmx::SelectionTreeElementPointer& sel,
 557                                       gmx_ana_index_t*                        g)
 558 {
 559     if (sel->child->evaluate)
 560     {
 561         sel->child->evaluate(data, sel->child, g);
 562     }
 563     sel->v.nr = sel->child->v.nr;
 564 }
 565
 566 /*!
 567  * \param[in] data Data for the current frame.
 568  * \param[in] sel  Selection element being evaluated.
 569  * \param[in] g    Group for which \p sel should be evaluated.
 570  * \returns   0 on success, a non-zero error code on error.
 571  *
 572  * If this is the first call for this frame, evaluates the child element
 573  * there should be exactly one in \p g.
 574  * The compiler has taken care that the child actually stores the evaluated
 575  * value in the value pointer of this element.
 576  * Assumes that \p g is persistent for the duration of the whole evaluation.
 577  *
 578  * This function is used as gmx::SelectionTreeElement::evaluate for
 579  * \ref SEL_SUBEXPR elements that have a static evaluation group, and hence do
 580  * not need full subexpression handling.
 581  */
 582 void _gmx_sel_evaluate_subexpr_staticeval(gmx_sel_evaluate_t*                     data,
 583                                           const gmx::SelectionTreeElementPointer& sel,
 584                                           gmx_ana_index_t*                        g)
 585 {
 586     if (sel->u.cgrp.isize == 0)
 587     {
 588         sel->child->evaluate(data, sel->child, g);
 589         sel->v.nr = sel->child->v.nr;
 590         if (!g)
 591         {
 592             sel->u.cgrp.isize = -1;
 593         }
 594         else
 595         {
 596             gmx_ana_index_set(&sel->u.cgrp, g->isize, g->index, 0);
 597         }
 598     }
 599 }
 600
 601 /*!
 602  * \param[in]  data  Data for the current frame.
 603  * \param[in]  sel   Selection element being evaluated.
 604  * \param[in]  g     Group for which \p sel should be evaluated.
 605  * \returns    0 on success, a non-zero error code on error.
 606  *
 607  * Finds the part of \p g for which the subexpression
 608  * has not yet been evaluated by comparing \p g to \p sel->u.cgrp.
 609  * If the part is not empty, the child expression is evaluated for this
 610  * part, and the results merged to the old values of the child.
 611  * The value of \p sel itself is undefined after the call.
 612  *
 613  * \todo
 614  * The call to gmx_ana_index_difference() can take quite a lot of unnecessary
 615  * time if the subexpression is evaluated either several times for the same
 616  * group or for completely distinct groups.
 617  * However, in the majority of cases, these situations occur when
 618  * _gmx_sel_evaluate_subexpr_staticeval() can be used, so this should not be a
 619  * major problem.
 620  */
 621 void _gmx_sel_evaluate_subexpr(gmx_sel_evaluate_t*                     data,
 622                                const gmx::SelectionTreeElementPointer& sel,
 623                                gmx_ana_index_t*                        g)
 624 {
 625     gmx_ana_index_t gmiss;
 626
 627     MempoolGroupReserver gmissreserver(data->mp);
 628     if (sel->u.cgrp.isize == 0)
 629     {
 630         {
 631             SelelemTemporaryValueAssigner assigner(sel->child, *sel);
 632             sel->child->evaluate(data, sel->child, g);
 633         }
 634         gmx_ana_index_copy(&sel->u.cgrp, g, false);
 635         gmiss.isize = 0;
 636     }
 637     else
 638     {
 639         gmissreserver.reserve(&gmiss, g->isize);
 640         gmx_ana_index_difference(&gmiss, g, &sel->u.cgrp);
 641     }
 642     if (gmiss.isize > 0)
 643     {
 644         MempoolSelelemReserver reserver(sel->child, gmiss.isize);
 645         /* Evaluate the missing values for the child */
 646         sel->child->evaluate(data, sel->child, &gmiss);
 647         /* Merge the missing values to the existing ones. */
 648         if (sel->v.type == GROUP_VALUE)
 649         {
 650             gmx_ana_index_merge(sel->v.u.g, sel->child->v.u.g, sel->v.u.g);
 651         }
 652         else
 653         {
 654             int i, j, k;
 655
 656             i = sel->u.cgrp.isize - 1;
 657             j = gmiss.isize - 1;
 658             /* TODO: This switch is kind of ugly, but it may be difficult to
 659              * do this portably without C++ templates. */
 660             switch (sel->v.type)
 661             {
 662                 case INT_VALUE:
 663                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 664                     {
 665                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 666                         {
 667                             sel->v.u.i[k] = sel->child->v.u.i[j--];
 668                         }
 669                         else
 670                         {
 671                             sel->v.u.i[k] = sel->v.u.i[i--];
 672                         }
 673                     }
 674                     break;
 675
 676                 case REAL_VALUE:
 677                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 678                     {
 679                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 680                         {
 681                             sel->v.u.r[k] = sel->child->v.u.r[j--];
 682                         }
 683                         else
 684                         {
 685                             sel->v.u.r[k] = sel->v.u.r[i--];
 686                         }
 687                     }
 688                     break;
 689
 690                 case STR_VALUE:
 691                     // Note: with the currently allowed syntax, this case is never
 692                     // reached.
 693                     for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
 694                     {
 695                         if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
 696                         {
 697                             sel->v.u.s[k] = sel->child->v.u.s[j--];
 698                         }
 699                         else
 700                         {
 701                             sel->v.u.s[k] = sel->v.u.s[i--];
 702                         }
 703                     }
 704                     break;
 705
 706                 case POS_VALUE:
 707                     /* TODO: Implement this */
 708                     GMX_THROW(gmx::NotImplementedError(
 709                             "position subexpressions not implemented properly"));
 710
 711                 case NO_VALUE:
 712                 case GROUP_VALUE: GMX_THROW(gmx::InternalError("Invalid subexpression type"));
 713             }
 714         }
 715         gmx_ana_index_merge(&sel->u.cgrp, &sel->u.cgrp, &gmiss);
 716     }
 717 }
 718
 719 /*!
 720  * \param[in] data Data for the current frame.
 721  * \param[in] sel Selection element being evaluated.
 722  * \param[in] g   Group for which \p sel should be evaluated.
 723  * \returns   0 for success.
 724  *
 725  * Sets the value pointers of the child and its child to point to the same
 726  * memory as the value pointer of this element to avoid copying, and then
 727  * evaluates evaluates the child.
 728  *
 729  * This function is used as gmx::SelectionTreeElement:evaluate for
 730  * \ref SEL_SUBEXPRREF elements for which the \ref SEL_SUBEXPR does not have
 731  * other references.
 732  */
 733 void _gmx_sel_evaluate_subexprref_simple(gmx_sel_evaluate_t*                     data,
 734                                          const gmx::SelectionTreeElementPointer& sel,
 735                                          gmx_ana_index_t*                        g)
 736 {
 737     if (g)
 738     {
 739         _gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
 740         _gmx_selvalue_setstore_alloc(&sel->child->child->v, sel->v.u.ptr, sel->child->child->v.nalloc);
 741         sel->child->evaluate(data, sel->child, g);
 742     }
 743     sel->v.nr = sel->child->v.nr;
 744     if (sel->u.param)
 745     {
 746         sel->u.param->val.nr = sel->v.nr;
 747         if (sel->u.param->nvalptr)
 748         {
 749             *sel->u.param->nvalptr = sel->u.param->val.nr;
 750         }
 751     }
 752 }
 753
 754 /*!
 755  * \param[in] data Data for the current frame.
 756  * \param[in] sel Selection element being evaluated.
 757  * \param[in] g   Group for which \p sel should be evaluated.
 758  * \returns   0 on success, a non-zero error code on error.
 759  *
 760  * If the value type is \ref POS_VALUE, the value of the child is simply
 761  * copied to set the value of \p sel (the child subexpression should
 762  * already have been evaluated by its root).
 763  * If the value type is something else, the child is evaluated for the
 764  * group \p g, and the value of the child is then copied.
 765  * There should be only one child element.
 766  *
 767  * This function is used as gmx::SelectionTreeElement::evaluate for
 768  * \ref SEL_SUBEXPRREF elements.
 769  */
 770 void _gmx_sel_evaluate_subexprref(gmx_sel_evaluate_t*                     data,
 771                                   const gmx::SelectionTreeElementPointer& sel,
 772                                   gmx_ana_index_t*                        g)
 773 {
 774     int i, j;
 775
 776     if (g != nullptr && sel->child->evaluate != nullptr)
 777     {
 778         sel->child->evaluate(data, sel->child, g);
 779     }
 780     const SelectionTreeElementPointer& expr = sel->child;
 781     switch (sel->v.type)
 782     {
 783         case INT_VALUE:
 784             if (!g)
 785             {
 786                 sel->v.nr = expr->v.nr;
 787                 memcpy(sel->v.u.i, expr->v.u.i, sel->v.nr * sizeof(*sel->v.u.i));
 788             }
 789             else
 790             {
 791                 sel->v.nr = g->isize;
 792                 /* Extract the values corresponding to g */
 793                 for (i = j = 0; i < g->isize; ++i, ++j)
 794                 {
 795                     while (sel->child->u.cgrp.index[j] < g->index[i])
 796                     {
 797                         ++j;
 798                     }
 799                     sel->v.u.i[i] = expr->v.u.i[j];
 800                 }
 801             }
 802             break;
 803
 804         case REAL_VALUE:
 805             if (!g)
 806             {
 807                 sel->v.nr = expr->v.nr;
 808                 memcpy(sel->v.u.r, expr->v.u.r, sel->v.nr * sizeof(*sel->v.u.r));
 809             }
 810             else
 811             {
 812                 sel->v.nr = g->isize;
 813                 /* Extract the values corresponding to g */
 814                 for (i = j = 0; i < g->isize; ++i, ++j)
 815                 {
 816                     while (sel->child->u.cgrp.index[j] < g->index[i])
 817                     {
 818                         ++j;
 819                     }
 820                     sel->v.u.r[i] = expr->v.u.r[j];
 821                 }
 822             }
 823             break;
 824
 825         case STR_VALUE:
 826             if (!g)
 827             {
 828                 sel->v.nr = expr->v.nr;
 829                 memcpy(sel->v.u.s, expr->v.u.s, sel->v.nr * sizeof(*sel->v.u.s));
 830             }
 831             else
 832             {
 833                 sel->v.nr = g->isize;
 834                 /* Extract the values corresponding to g */
 835                 for (i = j = 0; i < g->isize; ++i, ++j)
 836                 {
 837                     while (sel->child->u.cgrp.index[j] < g->index[i])
 838                     {
 839                         ++j;
 840                     }
 841                     sel->v.u.s[i] = expr->v.u.s[j];
 842                 }
 843             }
 844             break;
 845
 846         case POS_VALUE:
 847             /* Currently, there is no need to do anything fancy here,
 848              * but some future extensions may need a more flexible
 849              * implementation. */
 850             gmx_ana_pos_copy(sel->v.u.p, expr->v.u.p, false);
 851             break;
 852
 853         case GROUP_VALUE:
 854             if (!g)
 855             {
 856                 gmx_ana_index_copy(sel->v.u.g, expr->v.u.g, false);
 857             }
 858             else
 859             {
 860                 gmx_ana_index_intersection(sel->v.u.g, expr->v.u.g, g);
 861             }
 862             break;
 863
 864         default: /* should not be reached */
 865             GMX_THROW(gmx::InternalError("Invalid subexpression reference type"));
 866     }
 867     /* Store the number of values if needed */
 868     if (sel->u.param)
 869     {
 870         sel->u.param->val.nr = sel->v.nr;
 871         if (sel->u.param->nvalptr)
 872         {
 873             *sel->u.param->nvalptr = sel->u.param->val.nr;
 874         }
 875     }
 876 }
 877
 878 /********************************************************************
 879  * METHOD EXPRESSION EVALUATION
 880  ********************************************************************/
 881
 882 /*!
 883  * \param[in] data Data for the current frame.
 884  * \param[in] sel Selection element being evaluated.
 885  * \param[in] g   Group for which \p sel should be evaluated.
 886  * \returns   0 on success, a non-zero error code on error.
 887  *
 888  * Evaluates each child of a \ref SEL_EXPRESSION element.
 889  * The value of \p sel is not touched.
 890  *
 891  * This function is not used as gmx::SelectionTreeElement::evaluate,
 892  * but is used internally.
 893  */
 894 void _gmx_sel_evaluate_method_params(gmx_sel_evaluate_t*                     data,
 895                                      const gmx::SelectionTreeElementPointer& sel,
 896                                      gmx_ana_index_t*                        g)
 897 {
 898     SelectionTreeElementPointer child = sel->child;
 899     while (child)
 900     {
 901         if (child->evaluate && !(child->flags & SEL_EVALFRAME))
 902         {
 903             if (child->flags & SEL_ATOMVAL)
 904             {
 905                 child->evaluate(data, child, g);
 906             }
 907             else
 908             {
 909                 child->flags |= SEL_EVALFRAME;
 910                 child->evaluate(data, child, nullptr);
 911             }
 912         }
 913         child = child->next;
 914     }
 915 }
 916
 917 /*!
 918  * \param[in] data Data for the current frame.
 919  * \param[in] sel Selection element being evaluated.
 920  * \param[in] g   Group for which \p sel should be evaluated.
 921  * \returns   0 on success, a non-zero error code on error.
 922  *
 923  * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
 924  * to evaluate any parameter values.
 925  * If this is the first time this expression is evaluated for
 926  * the frame, sel_framefunc() callback is called if one is provided.
 927  * If a reference position calculation has been initialized for this element,
 928  * the positions are also updated, and sel_updatefunc_pos() is used to
 929  * evaluate the value. Otherwise, sel_updatefunc() is used.
 930  *
 931  * This function is used as gmx::SelectionTreeElement::evaluate for
 932  * \ref SEL_EXPRESSION elements.
 933  */
 934 void _gmx_sel_evaluate_method(gmx_sel_evaluate_t*                     data,
 935                               const gmx::SelectionTreeElementPointer& sel,
 936                               gmx_ana_index_t*                        g)
 937 {
 938     _gmx_sel_evaluate_method_params(data, sel, g);
 939     gmx::SelMethodEvalContext context(data->top, data->fr, data->pbc);
 940     if (sel->flags & SEL_INITFRAME)
 941     {
 942         sel->flags &= ~SEL_INITFRAME;
 943         sel->u.expr.method->init_frame(context, sel->u.expr.mdata);
 944     }
 945     if (sel->u.expr.pc)
 946     {
 947         gmx_ana_poscalc_update(sel->u.expr.pc, sel->u.expr.pos, g, data->fr, data->pbc);
 948         sel->u.expr.method->pupdate(context, sel->u.expr.pos, &sel->v, sel->u.expr.mdata);
 949         if ((sel->flags & SEL_ATOMVAL) && sel->v.nr < g->isize)
 950         {
 951             switch (sel->v.type)
 952             {
 953                 case REAL_VALUE:
 954                     expandValueForPositions(sel->v.u.r, &sel->v.nr, sel->u.expr.pos);
 955                     break;
 956                 default:
 957                     GMX_RELEASE_ASSERT(false, "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 }