src/gmxlib/selection/evaluate.c

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