src/gromacs/selection/evaluate.cpp

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