src/gromacs/selection/evaluate.cpp

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