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