src/gmxlib/selection/evaluate.c

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