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