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