2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements functions in evaluate.h.
40 * One of the major bottlenecks for selection performance is that all the
41 * evaluation is carried out for atoms.
42 * There are several cases when the evaluation could be done for residues
43 * or molecules instead, including keywords that select by residue and
44 * cases where residue centers are used as reference positions.
45 * Implementing this would require a mechanism for recognizing whether
46 * something can be evaluated by residue/molecule instead by atom, and
47 * converting selections by residue/molecule into selections by atom
50 * \author Teemu Murtola <teemu.murtola@gmail.com>
51 * \ingroup module_selection
61 #include "gromacs/math/utilities.h"
62 #include "gromacs/math/vec.h"
63 #include "gromacs/selection/indexutil.h"
64 #include "gromacs/selection/selection.h"
65 #include "gromacs/utility/exceptions.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/smalloc.h"
71 #include "selectioncollection_impl.h"
73 #include "selmethod.h"
75 using gmx::SelectionTreeElement;
76 using gmx::SelectionTreeElementPointer;
82 * Reserves memory for a selection element from the evaluation memory pool.
84 * This class implements RAII semantics for allocating memory for selection
85 * element values from a selection evaluation memory pool.
87 * \ingroup module_selection
89 class MempoolSelelemReserver
92 //! Constructs a reserver without initial reservation.
93 MempoolSelelemReserver() {}
95 * Constructs a reserver with initial reservation.
97 * \param[in,out] sel Selection element for which to reserve.
98 * \param[in] count Number of values to reserve.
102 MempoolSelelemReserver(const SelectionTreeElementPointer& sel, int count)
106 //! Frees any memory allocated using this reserver.
107 ~MempoolSelelemReserver()
111 sel_->mempoolRelease();
116 * Reserves memory for selection element values using this reserver.
118 * \param[in,out] sel Selection element for which to reserve.
119 * \param[in] count Number of values to reserve.
121 * Allocates space to store \p count output values in \p sel from the
122 * memory pool associated with \p sel, or from the heap if there is no
123 * memory pool. Type of values to allocate is automatically determined
126 void reserve(const SelectionTreeElementPointer& sel, int count)
128 GMX_RELEASE_ASSERT(!sel_, "Can only reserve one element with one instance");
129 sel->mempoolReserve(count);
134 SelectionTreeElementPointer sel_;
138 * Reserves memory for an index group from the evaluation memory pool.
140 * This class implements RAII semantics for allocating memory for an index
141 * group from a selection evaluation memory pool.
143 * \ingroup module_selection
145 class MempoolGroupReserver
149 * Creates a reserver associated with a given memory pool.
151 * \param mp Memory pool from which to reserve memory.
153 explicit MempoolGroupReserver(gmx_sel_mempool_t* mp) : mp_(mp), g_(nullptr) {}
154 //! Frees any memory allocated using this reserver.
155 ~MempoolGroupReserver()
159 _gmx_sel_mempool_free_group(mp_, g_);
164 * Reserves memory for an index group using this reserver.
166 * \param[in,out] g Index group to reserve.
167 * \param[in] count Number of atoms to reserve space for.
169 * Allocates memory from the memory pool to store \p count atoms in
172 void reserve(gmx_ana_index_t* g, int count)
174 GMX_RELEASE_ASSERT(g_ == nullptr, "Can only reserve one element with one instance");
175 _gmx_sel_mempool_alloc_group(mp_, g, count);
180 gmx_sel_mempool_t* mp_;
185 * Assigns a temporary value for a selection element.
187 * This class implements RAII semantics for temporarily assigning the value
188 * pointer of a selection element to point to a different location.
190 * \ingroup module_selection
192 class SelelemTemporaryValueAssigner
195 //! Constructs an assigner without an initial assignment.
196 SelelemTemporaryValueAssigner() : old_ptr_(nullptr), old_nalloc_(0) {}
198 * Constructs an assigner with an initial assignment.
200 * \param[in,out] sel Selection element for which to assign.
201 * \param[in] vsource Element to which \p sel values will point to.
205 SelelemTemporaryValueAssigner(const SelectionTreeElementPointer& sel, const SelectionTreeElement& vsource)
207 assign(sel, vsource);
209 //! Undoes any temporary assignment done using this assigner.
210 ~SelelemTemporaryValueAssigner()
214 _gmx_selvalue_setstore_alloc(&sel_->v, old_ptr_, old_nalloc_);
219 * Assigns a temporary value pointer.
221 * \param[in,out] sel Selection element for which to assign.
222 * \param[in] vsource Element to which \p sel values will point to.
224 * Assigns the value pointer in \p sel to point to the values in
225 * \p vsource, i.e., any access/modification to values in \p sel
226 * actually accesses values in \p vsource.
228 void assign(const SelectionTreeElementPointer& sel, const SelectionTreeElement& vsource)
230 GMX_RELEASE_ASSERT(!sel_, "Can only assign one element with one instance");
231 GMX_RELEASE_ASSERT(sel->v.type == vsource.v.type, "Mismatching selection value types");
232 _gmx_selvalue_getstore_and_release(&sel->v, &old_ptr_, &old_nalloc_);
233 _gmx_selvalue_setstore(&sel->v, vsource.v.u.ptr);
238 SelectionTreeElementPointer sel_;
244 * Expands a value array from one-per-position to one-per-atom.
246 * \param[in,out] value Array to expand.
247 * \param[in,out] nr Number of values in \p value.
248 * \param[in] pos Position data.
249 * \tparam T Value type of the array to expand.
251 * On input, \p value contains one value for each position in \p pos (and `*nr`
252 * must match). On output, \p value will contain a value for each atom used to
253 * evaluate `pos`: each input value is replicated to all atoms that make up the
254 * corresponding position.
255 * The operation is done in-place.
260 void expandValueForPositions(T value[], int* nr, gmx_ana_pos_t* pos)
262 GMX_RELEASE_ASSERT(*nr == pos->count(),
263 "Position update method did not return the correct number of values");
264 *nr = pos->m.mapb.nra;
265 // Loop over the positions in reverse order so that the expansion can be
266 // done in-place (assumes that each position has at least one atom, which
267 // should always be the case).
268 int outputIndex = pos->m.mapb.nra;
269 for (int i = pos->count() - 1; i >= 0; --i)
271 const int atomCount = pos->m.mapb.index[i + 1] - pos->m.mapb.index[i];
272 outputIndex -= atomCount;
273 GMX_ASSERT(outputIndex >= i, "In-place algorithm would overwrite data not yet used");
274 std::fill(&value[outputIndex], &value[outputIndex + atomCount], value[i]);
281 * \param[in] fp File handle to receive the output.
282 * \param[in] evalfunc Function pointer to print.
284 void _gmx_sel_print_evalfunc_name(FILE* fp, gmx::sel_evalfunc evalfunc)
290 else if (evalfunc == &_gmx_sel_evaluate_root)
294 else if (evalfunc == &_gmx_sel_evaluate_static)
296 fprintf(fp, "static");
298 else if (evalfunc == &_gmx_sel_evaluate_subexpr_simple)
300 fprintf(fp, "subexpr_simple");
302 else if (evalfunc == &_gmx_sel_evaluate_subexpr_staticeval)
304 fprintf(fp, "subexpr_staticeval");
306 else if (evalfunc == &_gmx_sel_evaluate_subexpr)
308 fprintf(fp, "subexpr");
310 else if (evalfunc == &_gmx_sel_evaluate_subexprref_simple)
312 fprintf(fp, "ref_simple");
314 else if (evalfunc == &_gmx_sel_evaluate_subexprref)
318 else if (evalfunc == &_gmx_sel_evaluate_method)
320 fprintf(fp, "method");
322 else if (evalfunc == &_gmx_sel_evaluate_modifier)
326 else if (evalfunc == &_gmx_sel_evaluate_not)
330 else if (evalfunc == &_gmx_sel_evaluate_and)
334 else if (evalfunc == &_gmx_sel_evaluate_or)
338 else if (evalfunc == &_gmx_sel_evaluate_arithmetic)
340 fprintf(fp, "arithmetic");
344 fprintf(fp, "%p", reinterpret_cast<void*>(evalfunc));
349 * \param[out] data Evaluation data structure to initialize.
350 * \param[in] mp Memory pool for intermediate evaluation values.
351 * \param[in] gall Index group with all the atoms.
352 * \param[in] top Topology structure for evaluation.
353 * \param[in] fr New frame for evaluation.
354 * \param[in] pbc New PBC information for evaluation.
356 void _gmx_sel_evaluate_init(gmx_sel_evaluate_t* data,
357 gmx_sel_mempool_t* mp,
358 gmx_ana_index_t* gall,
359 const gmx_mtop_t* top,
371 * Recursively initializes the flags for evaluation.
373 * \param[in,out] sel Selection element to clear.
375 * The \ref SEL_INITFRAME flag is set for \ref SEL_EXPRESSION elements whose
376 * method defines the \p init_frame callback (see sel_framefunc()), and
377 * cleared for other elements.
379 * The \ref SEL_EVALFRAME flag is cleared for all elements.
381 static void init_frame_eval(SelectionTreeElementPointer sel)
385 sel->flags &= ~(SEL_INITFRAME | SEL_EVALFRAME);
386 if (sel->type == SEL_EXPRESSION)
388 if (sel->u.expr.method && sel->u.expr.method->init_frame)
390 sel->flags |= SEL_INITFRAME;
393 if (sel->child && sel->type != SEL_SUBEXPRREF)
395 init_frame_eval(sel->child);
404 SelectionEvaluator::SelectionEvaluator() {}
407 * \param[in,out] coll The selection collection to evaluate.
408 * \param[in] fr Frame for which the evaluation should be carried out.
409 * \param[in] pbc PBC data, or NULL if no PBC should be used.
410 * \returns 0 on successful evaluation, a non-zero error code on error.
412 * This functions sets the global variables for topology, frame and PBC,
413 * clears some information in the selection to initialize the evaluation
414 * for a new frame, and evaluates \p sel and all the selections pointed by
415 * the \p next pointers of \p sel.
417 * This is the only function that user code should call if they want to
418 * evaluate a selection for a new frame.
420 void SelectionEvaluator::evaluate(SelectionCollection* coll, t_trxframe* fr, t_pbc* pbc)
422 gmx_ana_selcollection_t* sc = &coll->impl_->sc_;
423 gmx_sel_evaluate_t data;
425 _gmx_sel_evaluate_init(&data, sc->mempool, &sc->gall, sc->top, fr, pbc);
426 init_frame_eval(sc->root);
427 SelectionTreeElementPointer sel = sc->root;
430 /* Clear the evaluation group of subexpressions */
431 if (sel->child && sel->child->type == SEL_SUBEXPR && sel->child->evaluate != nullptr)
433 sel->child->u.cgrp.isize = 0;
434 /* Not strictly necessary, because the value will be overwritten
435 * during first evaluation of the subexpression anyways, but we
436 * clear the group for clarity. Note that this is _not_ done during
437 * compilation because of some additional complexities involved
438 * (see compiler.cpp), so it should not be relied upon in
439 * _gmx_sel_evaluate_subexpr(). */
440 if (sel->child->v.type == GROUP_VALUE)
442 sel->child->v.u.g->isize = 0;
447 sel->evaluate(&data, sel, nullptr);
451 /* Update selection information */
452 SelectionDataList::const_iterator isel;
453 for (isel = sc->sel.begin(); isel != sc->sel.end(); ++isel)
455 internal::SelectionData& sel = **isel;
456 sel.refreshMassesAndCharges(sc->top);
457 sel.updateCoveredFractionForFrame();
462 * \param[in,out] coll The selection collection to evaluate.
463 * \param[in] nframes Total number of frames.
465 void SelectionEvaluator::evaluateFinal(SelectionCollection* coll, int nframes)
467 gmx_ana_selcollection_t* sc = &coll->impl_->sc_;
469 SelectionDataList::const_iterator isel;
470 for (isel = sc->sel.begin(); isel != sc->sel.end(); ++isel)
472 internal::SelectionData& sel = **isel;
473 sel.restoreOriginalPositions(sc->top);
474 sel.computeAverageCoveredFraction(nframes);
481 * \param[in] data Data for the current frame.
482 * \param[in] sel Selection element being evaluated.
483 * \param[in] g Group for which \p sel should be evaluated.
484 * \returns 0 on success, a non-zero error code on error.
486 * Evaluates each child of \p sel in \p g.
488 void _gmx_sel_evaluate_children(gmx_sel_evaluate_t* data,
489 const gmx::SelectionTreeElementPointer& sel,
492 SelectionTreeElementPointer child = sel->child;
497 child->evaluate(data, child, g);
503 void _gmx_sel_evaluate_root(gmx_sel_evaluate_t* data,
504 const gmx::SelectionTreeElementPointer& sel,
505 gmx_ana_index_t* /* g */)
507 if (sel->u.cgrp.isize == 0 || !sel->child->evaluate)
512 sel->child->evaluate(data, sel->child, sel->u.cgrp.isize < 0 ? nullptr : &sel->u.cgrp);
515 void _gmx_sel_evaluate_static(gmx_sel_evaluate_t* /* data */,
516 const gmx::SelectionTreeElementPointer& sel,
519 if (sel->flags & SEL_UNSORTED)
521 // This only works if g contains all the atoms, but that is currently
522 // the only supported case.
523 gmx_ana_index_copy(sel->v.u.g, &sel->u.cgrp, false);
527 gmx_ana_index_intersection(sel->v.u.g, &sel->u.cgrp, g);
532 /*********************************************************************
533 * SUBEXPRESSION EVALUATION
534 *********************************************************************/
537 * \param[in] data Data for the current frame.
538 * \param[in] sel Selection element being evaluated.
539 * \param[in] g Group for which \p sel should be evaluated.
540 * \returns 0 on success, a non-zero error code on error.
542 * Evaluates the child element (there should be exactly one) in \p g.
543 * The compiler has taken care that the child actually stores the evaluated
544 * value in the value pointer of this element.
546 * This function is used as gmx::SelectionTreeElement::evaluate for
547 * \ref SEL_SUBEXPR elements that are used only once, and hence do not need
548 * full subexpression handling.
550 void _gmx_sel_evaluate_subexpr_simple(gmx_sel_evaluate_t* data,
551 const gmx::SelectionTreeElementPointer& sel,
554 if (sel->child->evaluate)
556 sel->child->evaluate(data, sel->child, g);
558 sel->v.nr = sel->child->v.nr;
562 * \param[in] data Data for the current frame.
563 * \param[in] sel Selection element being evaluated.
564 * \param[in] g Group for which \p sel should be evaluated.
565 * \returns 0 on success, a non-zero error code on error.
567 * If this is the first call for this frame, evaluates the child element
568 * there should be exactly one in \p g.
569 * The compiler has taken care that the child actually stores the evaluated
570 * value in the value pointer of this element.
571 * Assumes that \p g is persistent for the duration of the whole evaluation.
573 * This function is used as gmx::SelectionTreeElement::evaluate for
574 * \ref SEL_SUBEXPR elements that have a static evaluation group, and hence do
575 * not need full subexpression handling.
577 void _gmx_sel_evaluate_subexpr_staticeval(gmx_sel_evaluate_t* data,
578 const gmx::SelectionTreeElementPointer& sel,
581 if (sel->u.cgrp.isize == 0)
583 sel->child->evaluate(data, sel->child, g);
584 sel->v.nr = sel->child->v.nr;
587 sel->u.cgrp.isize = -1;
591 gmx_ana_index_set(&sel->u.cgrp, g->isize, g->index, 0);
597 * \param[in] data Data for the current frame.
598 * \param[in] sel Selection element being evaluated.
599 * \param[in] g Group for which \p sel should be evaluated.
600 * \returns 0 on success, a non-zero error code on error.
602 * Finds the part of \p g for which the subexpression
603 * has not yet been evaluated by comparing \p g to \p sel->u.cgrp.
604 * If the part is not empty, the child expression is evaluated for this
605 * part, and the results merged to the old values of the child.
606 * The value of \p sel itself is undefined after the call.
609 * The call to gmx_ana_index_difference() can take quite a lot of unnecessary
610 * time if the subexpression is evaluated either several times for the same
611 * group or for completely distinct groups.
612 * However, in the majority of cases, these situations occur when
613 * _gmx_sel_evaluate_subexpr_staticeval() can be used, so this should not be a
616 void _gmx_sel_evaluate_subexpr(gmx_sel_evaluate_t* data,
617 const gmx::SelectionTreeElementPointer& sel,
620 gmx_ana_index_t gmiss;
622 MempoolGroupReserver gmissreserver(data->mp);
623 if (sel->u.cgrp.isize == 0)
626 SelelemTemporaryValueAssigner assigner(sel->child, *sel);
627 sel->child->evaluate(data, sel->child, g);
629 gmx_ana_index_copy(&sel->u.cgrp, g, false);
634 gmissreserver.reserve(&gmiss, g->isize);
635 gmx_ana_index_difference(&gmiss, g, &sel->u.cgrp);
639 MempoolSelelemReserver reserver(sel->child, gmiss.isize);
640 /* Evaluate the missing values for the child */
641 sel->child->evaluate(data, sel->child, &gmiss);
642 /* Merge the missing values to the existing ones. */
643 if (sel->v.type == GROUP_VALUE)
645 gmx_ana_index_merge(sel->v.u.g, sel->child->v.u.g, sel->v.u.g);
651 i = sel->u.cgrp.isize - 1;
653 /* TODO: This switch is kind of ugly, but it may be difficult to
654 * do this portably without C++ templates. */
658 for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
660 if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
662 sel->v.u.i[k] = sel->child->v.u.i[j--];
666 sel->v.u.i[k] = sel->v.u.i[i--];
672 for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
674 if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
676 sel->v.u.r[k] = sel->child->v.u.r[j--];
680 sel->v.u.r[k] = sel->v.u.r[i--];
686 // Note: with the currently allowed syntax, this case is never
688 for (k = sel->u.cgrp.isize + gmiss.isize - 1; k >= 0; k--)
690 if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
692 sel->v.u.s[k] = sel->child->v.u.s[j--];
696 sel->v.u.s[k] = sel->v.u.s[i--];
702 /* TODO: Implement this */
703 GMX_THROW(gmx::NotImplementedError(
704 "position subexpressions not implemented properly"));
707 case GROUP_VALUE: GMX_THROW(gmx::InternalError("Invalid subexpression type"));
710 gmx_ana_index_merge(&sel->u.cgrp, &sel->u.cgrp, &gmiss);
715 * \param[in] data Data for the current frame.
716 * \param[in] sel Selection element being evaluated.
717 * \param[in] g Group for which \p sel should be evaluated.
718 * \returns 0 for success.
720 * Sets the value pointers of the child and its child to point to the same
721 * memory as the value pointer of this element to avoid copying, and then
722 * evaluates evaluates the child.
724 * This function is used as gmx::SelectionTreeElement:evaluate for
725 * \ref SEL_SUBEXPRREF elements for which the \ref SEL_SUBEXPR does not have
728 void _gmx_sel_evaluate_subexprref_simple(gmx_sel_evaluate_t* data,
729 const gmx::SelectionTreeElementPointer& sel,
734 _gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
735 _gmx_selvalue_setstore_alloc(&sel->child->child->v, sel->v.u.ptr, sel->child->child->v.nalloc);
736 sel->child->evaluate(data, sel->child, g);
738 sel->v.nr = sel->child->v.nr;
741 sel->u.param->val.nr = sel->v.nr;
742 if (sel->u.param->nvalptr)
744 *sel->u.param->nvalptr = sel->u.param->val.nr;
750 * \param[in] data Data for the current frame.
751 * \param[in] sel Selection element being evaluated.
752 * \param[in] g Group for which \p sel should be evaluated.
753 * \returns 0 on success, a non-zero error code on error.
755 * If the value type is \ref POS_VALUE, the value of the child is simply
756 * copied to set the value of \p sel (the child subexpression should
757 * already have been evaluated by its root).
758 * If the value type is something else, the child is evaluated for the
759 * group \p g, and the value of the child is then copied.
760 * There should be only one child element.
762 * This function is used as gmx::SelectionTreeElement::evaluate for
763 * \ref SEL_SUBEXPRREF elements.
765 void _gmx_sel_evaluate_subexprref(gmx_sel_evaluate_t* data,
766 const gmx::SelectionTreeElementPointer& sel,
771 if (g != nullptr && sel->child->evaluate != nullptr)
773 sel->child->evaluate(data, sel->child, g);
775 const SelectionTreeElementPointer& expr = sel->child;
781 sel->v.nr = expr->v.nr;
782 memcpy(sel->v.u.i, expr->v.u.i, sel->v.nr * sizeof(*sel->v.u.i));
786 sel->v.nr = g->isize;
787 /* Extract the values corresponding to g */
788 for (i = j = 0; i < g->isize; ++i, ++j)
790 while (sel->child->u.cgrp.index[j] < g->index[i])
794 sel->v.u.i[i] = expr->v.u.i[j];
802 sel->v.nr = expr->v.nr;
803 memcpy(sel->v.u.r, expr->v.u.r, sel->v.nr * sizeof(*sel->v.u.r));
807 sel->v.nr = g->isize;
808 /* Extract the values corresponding to g */
809 for (i = j = 0; i < g->isize; ++i, ++j)
811 while (sel->child->u.cgrp.index[j] < g->index[i])
815 sel->v.u.r[i] = expr->v.u.r[j];
823 sel->v.nr = expr->v.nr;
824 memcpy(sel->v.u.s, expr->v.u.s, sel->v.nr * sizeof(*sel->v.u.s));
828 sel->v.nr = g->isize;
829 /* Extract the values corresponding to g */
830 for (i = j = 0; i < g->isize; ++i, ++j)
832 while (sel->child->u.cgrp.index[j] < g->index[i])
836 sel->v.u.s[i] = expr->v.u.s[j];
842 /* Currently, there is no need to do anything fancy here,
843 * but some future extensions may need a more flexible
845 gmx_ana_pos_copy(sel->v.u.p, expr->v.u.p, false);
851 gmx_ana_index_copy(sel->v.u.g, expr->v.u.g, false);
855 gmx_ana_index_intersection(sel->v.u.g, expr->v.u.g, g);
859 default: /* should not be reached */
860 GMX_THROW(gmx::InternalError("Invalid subexpression reference type"));
862 /* Store the number of values if needed */
865 sel->u.param->val.nr = sel->v.nr;
866 if (sel->u.param->nvalptr)
868 *sel->u.param->nvalptr = sel->u.param->val.nr;
873 /********************************************************************
874 * METHOD EXPRESSION EVALUATION
875 ********************************************************************/
878 * \param[in] data Data for the current frame.
879 * \param[in] sel Selection element being evaluated.
880 * \param[in] g Group for which \p sel should be evaluated.
881 * \returns 0 on success, a non-zero error code on error.
883 * Evaluates each child of a \ref SEL_EXPRESSION element.
884 * The value of \p sel is not touched.
886 * This function is not used as gmx::SelectionTreeElement::evaluate,
887 * but is used internally.
889 void _gmx_sel_evaluate_method_params(gmx_sel_evaluate_t* data,
890 const gmx::SelectionTreeElementPointer& sel,
893 SelectionTreeElementPointer child = sel->child;
896 if (child->evaluate && !(child->flags & SEL_EVALFRAME))
898 if (child->flags & SEL_ATOMVAL)
900 child->evaluate(data, child, g);
904 child->flags |= SEL_EVALFRAME;
905 child->evaluate(data, child, nullptr);
913 * \param[in] data Data for the current frame.
914 * \param[in] sel Selection element being evaluated.
915 * \param[in] g Group for which \p sel should be evaluated.
916 * \returns 0 on success, a non-zero error code on error.
918 * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
919 * to evaluate any parameter values.
920 * If this is the first time this expression is evaluated for
921 * the frame, sel_framefunc() callback is called if one is provided.
922 * If a reference position calculation has been initialized for this element,
923 * the positions are also updated, and sel_updatefunc_pos() is used to
924 * evaluate the value. Otherwise, sel_updatefunc() is used.
926 * This function is used as gmx::SelectionTreeElement::evaluate for
927 * \ref SEL_EXPRESSION elements.
929 void _gmx_sel_evaluate_method(gmx_sel_evaluate_t* data,
930 const gmx::SelectionTreeElementPointer& sel,
933 _gmx_sel_evaluate_method_params(data, sel, g);
934 gmx::SelMethodEvalContext context(data->top, data->fr, data->pbc);
935 if (sel->flags & SEL_INITFRAME)
937 sel->flags &= ~SEL_INITFRAME;
938 sel->u.expr.method->init_frame(context, sel->u.expr.mdata);
942 gmx_ana_poscalc_update(sel->u.expr.pc, sel->u.expr.pos, g, data->fr, data->pbc);
943 sel->u.expr.method->pupdate(context, sel->u.expr.pos, &sel->v, sel->u.expr.mdata);
944 if ((sel->flags & SEL_ATOMVAL) && sel->v.nr < g->isize)
949 expandValueForPositions(sel->v.u.r, &sel->v.nr, sel->u.expr.pos);
952 GMX_RELEASE_ASSERT(false,
953 "Unimplemented value type for position update method");
959 sel->u.expr.method->update(context, g, &sel->v, sel->u.expr.mdata);
964 * \param[in] data Data for the current frame.
965 * \param[in] sel Selection element being evaluated.
966 * \param[in] g Group for which \p sel should be evaluated.
967 * \returns 0 on success, a non-zero error code on error.
969 * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
970 * to evaluate any parameter values.
971 * If this is the first time this expression is evaluated for
972 * the frame, sel_framefunc() callback is called if one is provided.
973 * The modifier is then evaluated using sel_updatefunc_pos().
975 * This function is used as gmx::SelectionTreeElement::evaluate for
976 * \ref SEL_MODIFIER elements.
978 void _gmx_sel_evaluate_modifier(gmx_sel_evaluate_t* data,
979 const gmx::SelectionTreeElementPointer& sel,
982 _gmx_sel_evaluate_method_params(data, sel, g);
983 gmx::SelMethodEvalContext context(data->top, data->fr, data->pbc);
984 if (sel->flags & SEL_INITFRAME)
986 sel->flags &= ~SEL_INITFRAME;
987 sel->u.expr.method->init_frame(context, sel->u.expr.mdata);
989 if (sel->child && sel->child->v.type != POS_VALUE)
991 GMX_THROW(gmx::NotImplementedError("Non-position valued modifiers not implemented"));
993 sel->u.expr.method->pupdate(context, nullptr, &sel->v, sel->u.expr.mdata);
997 /********************************************************************
998 * BOOLEAN EXPRESSION EVALUATION
999 ********************************************************************/
1002 * \param[in] data Data for the current frame.
1003 * \param[in] sel Selection element being evaluated.
1004 * \param[in] g Group for which \p sel should be evaluated.
1005 * \returns 0 on success, a non-zero error code on error.
1007 * Evaluates the child element (there should be only one) in the group
1008 * \p g, and then sets the value of \p sel to the complement of the
1011 * This function is used as gmx::SelectionTreeElement::evaluate for
1012 * \ref SEL_BOOLEAN elements with \ref BOOL_NOT.
1014 void _gmx_sel_evaluate_not(gmx_sel_evaluate_t* data,
1015 const gmx::SelectionTreeElementPointer& sel,
1018 MempoolSelelemReserver reserver(sel->child, g->isize);
1019 sel->child->evaluate(data, sel->child, g);
1020 gmx_ana_index_difference(sel->v.u.g, g, sel->child->v.u.g);
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.
1029 * Short-circuiting evaluation of logical AND expressions.
1031 * Starts by evaluating the first child element in the group \p g.
1032 * The each following child element is evaluated in the intersection
1033 * of all the previous values until all children have been evaluated
1034 * or the intersection becomes empty.
1035 * The value of \p sel is set to the intersection of all the (evaluated)
1038 * If the first child does not have an evaluation function, it is skipped
1039 * and the evaluation is started at the second child.
1040 * This happens if the first child is a constant expression and during
1041 * compilation it was detected that the evaluation group is always a subset
1042 * of the constant group
1043 * (currently, the compiler never detects this).
1045 * This function is used as gmx::SelectionTreeElement::evaluate for
1046 * \ref SEL_BOOLEAN elements with \ref BOOL_AND.
1048 void _gmx_sel_evaluate_and(gmx_sel_evaluate_t* data,
1049 const gmx::SelectionTreeElementPointer& sel,
1052 SelectionTreeElementPointer child = sel->child;
1053 /* Skip the first child if it does not have an evaluation function. */
1054 if (!child->evaluate)
1056 child = child->next;
1058 /* Evaluate the first child */
1060 MempoolSelelemReserver reserver(child, g->isize);
1061 child->evaluate(data, child, g);
1062 gmx_ana_index_copy(sel->v.u.g, child->v.u.g, false);
1064 child = child->next;
1065 while (child && sel->v.u.g->isize > 0)
1067 MempoolSelelemReserver reserver(child, sel->v.u.g->isize);
1068 child->evaluate(data, child, sel->v.u.g);
1069 gmx_ana_index_intersection(sel->v.u.g, sel->v.u.g, child->v.u.g);
1070 child = child->next;
1075 * \param[in] data Data for the current frame.
1076 * \param[in] sel Selection element being evaluated.
1077 * \param[in] g Group for which \p sel should be evaluated.
1078 * \returns 0 on success, a non-zero error code on error.
1080 * Short-circuiting evaluation of logical OR expressions.
1082 * Starts by evaluating the first child element in the group \p g.
1083 * For each subsequent child, finds the part of \p g that is not
1084 * included the value of any previous child, and evaluates the child
1085 * in that group until the last child is evaluated or all of \p g
1086 * is included in some child value.
1087 * The value of \p sel is set to the union of all the (evaluated)
1090 * If the first child does not have an evaluation function, its value is
1091 * used without evaluation.
1092 * This happens if the first child is a constant expression, the selection
1093 * has been compiled, and the evaluation group is the same for each frame.
1094 * In this case, the compiler has taken care of that the child value is a
1095 * subset of \p g, making it unnecessary to evaluate it.
1097 * This function is used as gmx::SelectionTreeElement::evaluate for
1098 * \ref SEL_BOOLEAN elements with \ref BOOL_OR.
1100 void _gmx_sel_evaluate_or(gmx_sel_evaluate_t* data,
1101 const gmx::SelectionTreeElementPointer& sel,
1104 gmx_ana_index_t tmp, tmp2;
1106 SelectionTreeElementPointer child = sel->child;
1107 if (child->evaluate)
1109 MempoolSelelemReserver reserver(child, g->isize);
1110 child->evaluate(data, child, g);
1111 gmx_ana_index_partition(sel->v.u.g, &tmp, g, child->v.u.g);
1115 gmx_ana_index_partition(sel->v.u.g, &tmp, g, child->v.u.g);
1117 child = child->next;
1118 while (child && tmp.isize > 0)
1121 MempoolSelelemReserver reserver(child, tmp.isize);
1122 child->evaluate(data, child, &tmp);
1123 gmx_ana_index_partition(&tmp, &tmp2, &tmp, child->v.u.g);
1125 sel->v.u.g->isize += tmp.isize;
1126 tmp.isize = tmp2.isize;
1127 tmp.index = tmp2.index;
1128 child = child->next;
1130 gmx_ana_index_sort(sel->v.u.g);
1134 /********************************************************************
1135 * ARITHMETIC EVALUATION
1136 ********************************************************************/
1139 * \param[in] data Data for the current frame.
1140 * \param[in] sel Selection element being evaluated.
1141 * \param[in] g Group for which \p sel should be evaluated.
1142 * \returns 0 on success, a non-zero error code on error.
1144 void _gmx_sel_evaluate_arithmetic(gmx_sel_evaluate_t* data,
1145 const gmx::SelectionTreeElementPointer& sel,
1149 real lval, rval = 0., val = 0.;
1151 const SelectionTreeElementPointer& left = sel->child;
1152 const SelectionTreeElementPointer& right = left->next;
1154 SelelemTemporaryValueAssigner assigner;
1155 MempoolSelelemReserver reserver;
1158 assigner.assign(left, *sel);
1161 reserver.reserve(right, g->isize);
1164 else if (right && right->mempool)
1166 assigner.assign(right, *sel);
1168 _gmx_sel_evaluate_children(data, sel, g);
1170 n = (sel->flags & SEL_SINGLEVAL) ? 1 : g->isize;
1173 bool bArithNeg = (sel->u.arith.type == ARITH_NEG);
1174 GMX_ASSERT(right || bArithNeg, "Right operand cannot be null except for negations");
1175 for (i = i1 = i2 = 0; i < n; ++i)
1177 lval = left->v.u.r[i1];
1180 rval = right->v.u.r[i2];
1182 switch (sel->u.arith.type)
1184 case ARITH_PLUS: val = lval + rval; break;
1185 case ARITH_MINUS: val = lval - rval; break;
1186 case ARITH_NEG: val = -lval; break;
1187 case ARITH_MULT: val = lval * rval; break;
1188 case ARITH_DIV: val = lval / rval; break;
1189 case ARITH_EXP: val = pow(lval, rval); break;
1191 sel->v.u.r[i] = val;
1192 if (!(left->flags & SEL_SINGLEVAL))
1196 if (!bArithNeg && !(right->flags & SEL_SINGLEVAL))