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