2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team.
5 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
6 * Copyright (c) 2019,2020, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 * Implements functions in evaluate.h.
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
52 * \author Teemu Murtola <teemu.murtola@gmail.com>
53 * \ingroup module_selection
63 #include "gromacs/math/utilities.h"
64 #include "gromacs/math/vec.h"
65 #include "gromacs/selection/indexutil.h"
66 #include "gromacs/selection/selection.h"
67 #include "gromacs/utility/exceptions.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/smalloc.h"
73 #include "selectioncollection_impl.h"
75 #include "selmethod.h"
77 using gmx::SelectionTreeElement
;
78 using gmx::SelectionTreeElementPointer
;
84 * Reserves memory for a selection element from the evaluation memory pool.
86 * This class implements RAII semantics for allocating memory for selection
87 * element values from a selection evaluation memory pool.
89 * \ingroup module_selection
91 class MempoolSelelemReserver
94 //! Constructs a reserver without initial reservation.
95 MempoolSelelemReserver() {}
97 * Constructs a reserver with initial reservation.
99 * \param[in,out] sel Selection element for which to reserve.
100 * \param[in] count Number of values to reserve.
104 MempoolSelelemReserver(const SelectionTreeElementPointer
& sel
, int count
)
108 //! Frees any memory allocated using this reserver.
109 ~MempoolSelelemReserver()
113 sel_
->mempoolRelease();
118 * Reserves memory for selection element values using this reserver.
120 * \param[in,out] sel Selection element for which to reserve.
121 * \param[in] count Number of values to reserve.
123 * Allocates space to store \p count output values in \p sel from the
124 * memory pool associated with \p sel, or from the heap if there is no
125 * memory pool. Type of values to allocate is automatically determined
128 void reserve(const SelectionTreeElementPointer
& sel
, int count
)
130 GMX_RELEASE_ASSERT(!sel_
, "Can only reserve one element with one instance");
131 sel
->mempoolReserve(count
);
136 SelectionTreeElementPointer sel_
;
140 * Reserves memory for an index group from the evaluation memory pool.
142 * This class implements RAII semantics for allocating memory for an index
143 * group from a selection evaluation memory pool.
145 * \ingroup module_selection
147 class MempoolGroupReserver
151 * Creates a reserver associated with a given memory pool.
153 * \param mp Memory pool from which to reserve memory.
155 explicit MempoolGroupReserver(gmx_sel_mempool_t
* mp
) : mp_(mp
), g_(nullptr) {}
156 //! Frees any memory allocated using this reserver.
157 ~MempoolGroupReserver()
161 _gmx_sel_mempool_free_group(mp_
, g_
);
166 * Reserves memory for an index group using this reserver.
168 * \param[in,out] g Index group to reserve.
169 * \param[in] count Number of atoms to reserve space for.
171 * Allocates memory from the memory pool to store \p count atoms in
174 void reserve(gmx_ana_index_t
* g
, int count
)
176 GMX_RELEASE_ASSERT(g_
== nullptr, "Can only reserve one element with one instance");
177 _gmx_sel_mempool_alloc_group(mp_
, g
, count
);
182 gmx_sel_mempool_t
* mp_
;
187 * Assigns a temporary value for a selection element.
189 * This class implements RAII semantics for temporarily assigning the value
190 * pointer of a selection element to point to a different location.
192 * \ingroup module_selection
194 class SelelemTemporaryValueAssigner
197 //! Constructs an assigner without an initial assignment.
198 SelelemTemporaryValueAssigner() : old_ptr_(nullptr), old_nalloc_(0) {}
200 * Constructs an assigner with an initial assignment.
202 * \param[in,out] sel Selection element for which to assign.
203 * \param[in] vsource Element to which \p sel values will point to.
207 SelelemTemporaryValueAssigner(const SelectionTreeElementPointer
& sel
, const SelectionTreeElement
& vsource
)
209 assign(sel
, vsource
);
211 //! Undoes any temporary assignment done using this assigner.
212 ~SelelemTemporaryValueAssigner()
216 _gmx_selvalue_setstore_alloc(&sel_
->v
, old_ptr_
, old_nalloc_
);
221 * Assigns a temporary value pointer.
223 * \param[in,out] sel Selection element for which to assign.
224 * \param[in] vsource Element to which \p sel values will point to.
226 * Assigns the value pointer in \p sel to point to the values in
227 * \p vsource, i.e., any access/modification to values in \p sel
228 * actually accesses values in \p vsource.
230 void assign(const SelectionTreeElementPointer
& sel
, const SelectionTreeElement
& vsource
)
232 GMX_RELEASE_ASSERT(!sel_
, "Can only assign one element with one instance");
233 GMX_RELEASE_ASSERT(sel
->v
.type
== vsource
.v
.type
, "Mismatching selection value types");
234 _gmx_selvalue_getstore_and_release(&sel
->v
, &old_ptr_
, &old_nalloc_
);
235 _gmx_selvalue_setstore(&sel
->v
, vsource
.v
.u
.ptr
);
240 SelectionTreeElementPointer sel_
;
246 * Expands a value array from one-per-position to one-per-atom.
248 * \param[in,out] value Array to expand.
249 * \param[in,out] nr Number of values in \p value.
250 * \param[in] pos Position data.
251 * \tparam T Value type of the array to expand.
253 * On input, \p value contains one value for each position in \p pos (and `*nr`
254 * must match). On output, \p value will contain a value for each atom used to
255 * evaluate `pos`: each input value is replicated to all atoms that make up the
256 * corresponding position.
257 * The operation is done in-place.
262 void expandValueForPositions(T value
[], int* nr
, gmx_ana_pos_t
* pos
)
264 GMX_RELEASE_ASSERT(*nr
== pos
->count(),
265 "Position update method did not return the correct number of values");
266 *nr
= pos
->m
.mapb
.nra
;
267 // Loop over the positions in reverse order so that the expansion can be
268 // done in-place (assumes that each position has at least one atom, which
269 // should always be the case).
270 int outputIndex
= pos
->m
.mapb
.nra
;
271 for (int i
= pos
->count() - 1; i
>= 0; --i
)
273 const int atomCount
= pos
->m
.mapb
.index
[i
+ 1] - pos
->m
.mapb
.index
[i
];
274 outputIndex
-= atomCount
;
275 GMX_ASSERT(outputIndex
>= i
, "In-place algorithm would overwrite data not yet used");
276 std::fill(&value
[outputIndex
], &value
[outputIndex
+ atomCount
], value
[i
]);
283 * \param[in] fp File handle to receive the output.
284 * \param[in] evalfunc Function pointer to print.
286 void _gmx_sel_print_evalfunc_name(FILE* fp
, gmx::sel_evalfunc evalfunc
)
292 else if (evalfunc
== &_gmx_sel_evaluate_root
)
296 else if (evalfunc
== &_gmx_sel_evaluate_static
)
298 fprintf(fp
, "static");
300 else if (evalfunc
== &_gmx_sel_evaluate_subexpr_simple
)
302 fprintf(fp
, "subexpr_simple");
304 else if (evalfunc
== &_gmx_sel_evaluate_subexpr_staticeval
)
306 fprintf(fp
, "subexpr_staticeval");
308 else if (evalfunc
== &_gmx_sel_evaluate_subexpr
)
310 fprintf(fp
, "subexpr");
312 else if (evalfunc
== &_gmx_sel_evaluate_subexprref_simple
)
314 fprintf(fp
, "ref_simple");
316 else if (evalfunc
== &_gmx_sel_evaluate_subexprref
)
320 else if (evalfunc
== &_gmx_sel_evaluate_method
)
322 fprintf(fp
, "method");
324 else if (evalfunc
== &_gmx_sel_evaluate_modifier
)
328 else if (evalfunc
== &_gmx_sel_evaluate_not
)
332 else if (evalfunc
== &_gmx_sel_evaluate_and
)
336 else if (evalfunc
== &_gmx_sel_evaluate_or
)
340 else if (evalfunc
== &_gmx_sel_evaluate_arithmetic
)
342 fprintf(fp
, "arithmetic");
346 fprintf(fp
, "%p", reinterpret_cast<void*>(evalfunc
));
351 * \param[out] data Evaluation data structure to initialize.
352 * \param[in] mp Memory pool for intermediate evaluation values.
353 * \param[in] gall Index group with all the atoms.
354 * \param[in] top Topology structure for evaluation.
355 * \param[in] fr New frame for evaluation.
356 * \param[in] pbc New PBC information for evaluation.
358 void _gmx_sel_evaluate_init(gmx_sel_evaluate_t
* data
,
359 gmx_sel_mempool_t
* mp
,
360 gmx_ana_index_t
* gall
,
361 const gmx_mtop_t
* top
,
373 * Recursively initializes the flags for evaluation.
375 * \param[in,out] sel Selection element to clear.
377 * The \ref SEL_INITFRAME flag is set for \ref SEL_EXPRESSION elements whose
378 * method defines the \p init_frame callback (see sel_framefunc()), and
379 * cleared for other elements.
381 * The \ref SEL_EVALFRAME flag is cleared for all elements.
383 static void init_frame_eval(SelectionTreeElementPointer sel
)
387 sel
->flags
&= ~(SEL_INITFRAME
| SEL_EVALFRAME
);
388 if (sel
->type
== SEL_EXPRESSION
)
390 if (sel
->u
.expr
.method
&& sel
->u
.expr
.method
->init_frame
)
392 sel
->flags
|= SEL_INITFRAME
;
395 if (sel
->child
&& sel
->type
!= SEL_SUBEXPRREF
)
397 init_frame_eval(sel
->child
);
406 SelectionEvaluator::SelectionEvaluator() {}
409 * \param[in,out] coll The selection collection to evaluate.
410 * \param[in] fr Frame for which the evaluation should be carried out.
411 * \param[in] pbc PBC data, or NULL if no PBC should be used.
412 * \returns 0 on successful evaluation, a non-zero error code on error.
414 * This functions sets the global variables for topology, frame and PBC,
415 * clears some information in the selection to initialize the evaluation
416 * for a new frame, and evaluates \p sel and all the selections pointed by
417 * the \p next pointers of \p sel.
419 * This is the only function that user code should call if they want to
420 * evaluate a selection for a new frame.
422 void SelectionEvaluator::evaluate(SelectionCollection
* coll
, t_trxframe
* fr
, t_pbc
* pbc
)
424 gmx_ana_selcollection_t
* sc
= &coll
->impl_
->sc_
;
425 gmx_sel_evaluate_t data
;
427 _gmx_sel_evaluate_init(&data
, sc
->mempool
, &sc
->gall
, sc
->top
, fr
, pbc
);
428 init_frame_eval(sc
->root
);
429 SelectionTreeElementPointer sel
= sc
->root
;
432 /* Clear the evaluation group of subexpressions */
433 if (sel
->child
&& sel
->child
->type
== SEL_SUBEXPR
&& sel
->child
->evaluate
!= nullptr)
435 sel
->child
->u
.cgrp
.isize
= 0;
436 /* Not strictly necessary, because the value will be overwritten
437 * during first evaluation of the subexpression anyways, but we
438 * clear the group for clarity. Note that this is _not_ done during
439 * compilation because of some additional complexities involved
440 * (see compiler.cpp), so it should not be relied upon in
441 * _gmx_sel_evaluate_subexpr(). */
442 if (sel
->child
->v
.type
== GROUP_VALUE
)
444 sel
->child
->v
.u
.g
->isize
= 0;
449 sel
->evaluate(&data
, sel
, nullptr);
453 /* Update selection information */
454 SelectionDataList::const_iterator isel
;
455 for (isel
= sc
->sel
.begin(); isel
!= sc
->sel
.end(); ++isel
)
457 internal::SelectionData
& sel
= **isel
;
458 sel
.refreshMassesAndCharges(sc
->top
);
459 sel
.updateCoveredFractionForFrame();
464 * \param[in,out] coll The selection collection to evaluate.
465 * \param[in] nframes Total number of frames.
467 void SelectionEvaluator::evaluateFinal(SelectionCollection
* coll
, int nframes
)
469 gmx_ana_selcollection_t
* sc
= &coll
->impl_
->sc_
;
471 SelectionDataList::const_iterator isel
;
472 for (isel
= sc
->sel
.begin(); isel
!= sc
->sel
.end(); ++isel
)
474 internal::SelectionData
& sel
= **isel
;
475 sel
.restoreOriginalPositions(sc
->top
);
476 sel
.computeAverageCoveredFraction(nframes
);
483 * \param[in] data Data for the current frame.
484 * \param[in] sel Selection element being evaluated.
485 * \param[in] g Group for which \p sel should be evaluated.
486 * \returns 0 on success, a non-zero error code on error.
488 * Evaluates each child of \p sel in \p g.
490 void _gmx_sel_evaluate_children(gmx_sel_evaluate_t
* data
,
491 const gmx::SelectionTreeElementPointer
& sel
,
494 SelectionTreeElementPointer child
= sel
->child
;
499 child
->evaluate(data
, child
, g
);
505 void _gmx_sel_evaluate_root(gmx_sel_evaluate_t
* data
,
506 const gmx::SelectionTreeElementPointer
& sel
,
507 gmx_ana_index_t
* /* g */)
509 if (sel
->u
.cgrp
.isize
== 0 || !sel
->child
->evaluate
)
514 sel
->child
->evaluate(data
, sel
->child
, sel
->u
.cgrp
.isize
< 0 ? nullptr : &sel
->u
.cgrp
);
517 void _gmx_sel_evaluate_static(gmx_sel_evaluate_t
* /* data */,
518 const gmx::SelectionTreeElementPointer
& sel
,
521 if (sel
->flags
& SEL_UNSORTED
)
523 // This only works if g contains all the atoms, but that is currently
524 // the only supported case.
525 gmx_ana_index_copy(sel
->v
.u
.g
, &sel
->u
.cgrp
, false);
529 gmx_ana_index_intersection(sel
->v
.u
.g
, &sel
->u
.cgrp
, g
);
534 /*********************************************************************
535 * SUBEXPRESSION EVALUATION
536 *********************************************************************/
539 * \param[in] data Data for the current frame.
540 * \param[in] sel Selection element being evaluated.
541 * \param[in] g Group for which \p sel should be evaluated.
542 * \returns 0 on success, a non-zero error code on error.
544 * Evaluates the child element (there should be exactly one) in \p g.
545 * The compiler has taken care that the child actually stores the evaluated
546 * value in the value pointer of this element.
548 * This function is used as gmx::SelectionTreeElement::evaluate for
549 * \ref SEL_SUBEXPR elements that are used only once, and hence do not need
550 * full subexpression handling.
552 void _gmx_sel_evaluate_subexpr_simple(gmx_sel_evaluate_t
* data
,
553 const gmx::SelectionTreeElementPointer
& sel
,
556 if (sel
->child
->evaluate
)
558 sel
->child
->evaluate(data
, sel
->child
, g
);
560 sel
->v
.nr
= sel
->child
->v
.nr
;
564 * \param[in] data Data for the current frame.
565 * \param[in] sel Selection element being evaluated.
566 * \param[in] g Group for which \p sel should be evaluated.
567 * \returns 0 on success, a non-zero error code on error.
569 * If this is the first call for this frame, evaluates the child element
570 * there should be exactly one in \p g.
571 * The compiler has taken care that the child actually stores the evaluated
572 * value in the value pointer of this element.
573 * Assumes that \p g is persistent for the duration of the whole evaluation.
575 * This function is used as gmx::SelectionTreeElement::evaluate for
576 * \ref SEL_SUBEXPR elements that have a static evaluation group, and hence do
577 * not need full subexpression handling.
579 void _gmx_sel_evaluate_subexpr_staticeval(gmx_sel_evaluate_t
* data
,
580 const gmx::SelectionTreeElementPointer
& sel
,
583 if (sel
->u
.cgrp
.isize
== 0)
585 sel
->child
->evaluate(data
, sel
->child
, g
);
586 sel
->v
.nr
= sel
->child
->v
.nr
;
589 sel
->u
.cgrp
.isize
= -1;
593 gmx_ana_index_set(&sel
->u
.cgrp
, g
->isize
, g
->index
, 0);
599 * \param[in] data Data for the current frame.
600 * \param[in] sel Selection element being evaluated.
601 * \param[in] g Group for which \p sel should be evaluated.
602 * \returns 0 on success, a non-zero error code on error.
604 * Finds the part of \p g for which the subexpression
605 * has not yet been evaluated by comparing \p g to \p sel->u.cgrp.
606 * If the part is not empty, the child expression is evaluated for this
607 * part, and the results merged to the old values of the child.
608 * The value of \p sel itself is undefined after the call.
611 * The call to gmx_ana_index_difference() can take quite a lot of unnecessary
612 * time if the subexpression is evaluated either several times for the same
613 * group or for completely distinct groups.
614 * However, in the majority of cases, these situations occur when
615 * _gmx_sel_evaluate_subexpr_staticeval() can be used, so this should not be a
618 void _gmx_sel_evaluate_subexpr(gmx_sel_evaluate_t
* data
,
619 const gmx::SelectionTreeElementPointer
& sel
,
622 gmx_ana_index_t gmiss
;
624 MempoolGroupReserver
gmissreserver(data
->mp
);
625 if (sel
->u
.cgrp
.isize
== 0)
628 SelelemTemporaryValueAssigner
assigner(sel
->child
, *sel
);
629 sel
->child
->evaluate(data
, sel
->child
, g
);
631 gmx_ana_index_copy(&sel
->u
.cgrp
, g
, false);
636 gmissreserver
.reserve(&gmiss
, g
->isize
);
637 gmx_ana_index_difference(&gmiss
, g
, &sel
->u
.cgrp
);
641 MempoolSelelemReserver
reserver(sel
->child
, gmiss
.isize
);
642 /* Evaluate the missing values for the child */
643 sel
->child
->evaluate(data
, sel
->child
, &gmiss
);
644 /* Merge the missing values to the existing ones. */
645 if (sel
->v
.type
== GROUP_VALUE
)
647 gmx_ana_index_merge(sel
->v
.u
.g
, sel
->child
->v
.u
.g
, sel
->v
.u
.g
);
653 i
= sel
->u
.cgrp
.isize
- 1;
655 /* TODO: This switch is kind of ugly, but it may be difficult to
656 * do this portably without C++ templates. */
660 for (k
= sel
->u
.cgrp
.isize
+ gmiss
.isize
- 1; k
>= 0; k
--)
662 if (i
< 0 || (j
>= 0 && sel
->u
.cgrp
.index
[i
] < gmiss
.index
[j
]))
664 sel
->v
.u
.i
[k
] = sel
->child
->v
.u
.i
[j
--];
668 sel
->v
.u
.i
[k
] = sel
->v
.u
.i
[i
--];
674 for (k
= sel
->u
.cgrp
.isize
+ gmiss
.isize
- 1; k
>= 0; k
--)
676 if (i
< 0 || (j
>= 0 && sel
->u
.cgrp
.index
[i
] < gmiss
.index
[j
]))
678 sel
->v
.u
.r
[k
] = sel
->child
->v
.u
.r
[j
--];
682 sel
->v
.u
.r
[k
] = sel
->v
.u
.r
[i
--];
688 // Note: with the currently allowed syntax, this case is never
690 for (k
= sel
->u
.cgrp
.isize
+ gmiss
.isize
- 1; k
>= 0; k
--)
692 if (i
< 0 || (j
>= 0 && sel
->u
.cgrp
.index
[i
] < gmiss
.index
[j
]))
694 sel
->v
.u
.s
[k
] = sel
->child
->v
.u
.s
[j
--];
698 sel
->v
.u
.s
[k
] = sel
->v
.u
.s
[i
--];
704 /* TODO: Implement this */
705 GMX_THROW(gmx::NotImplementedError(
706 "position subexpressions not implemented properly"));
709 case GROUP_VALUE
: GMX_THROW(gmx::InternalError("Invalid subexpression type"));
712 gmx_ana_index_merge(&sel
->u
.cgrp
, &sel
->u
.cgrp
, &gmiss
);
717 * \param[in] data Data for the current frame.
718 * \param[in] sel Selection element being evaluated.
719 * \param[in] g Group for which \p sel should be evaluated.
720 * \returns 0 for success.
722 * Sets the value pointers of the child and its child to point to the same
723 * memory as the value pointer of this element to avoid copying, and then
724 * evaluates evaluates the child.
726 * This function is used as gmx::SelectionTreeElement:evaluate for
727 * \ref SEL_SUBEXPRREF elements for which the \ref SEL_SUBEXPR does not have
730 void _gmx_sel_evaluate_subexprref_simple(gmx_sel_evaluate_t
* data
,
731 const gmx::SelectionTreeElementPointer
& sel
,
736 _gmx_selvalue_setstore(&sel
->child
->v
, sel
->v
.u
.ptr
);
737 _gmx_selvalue_setstore_alloc(&sel
->child
->child
->v
, sel
->v
.u
.ptr
, sel
->child
->child
->v
.nalloc
);
738 sel
->child
->evaluate(data
, sel
->child
, g
);
740 sel
->v
.nr
= sel
->child
->v
.nr
;
743 sel
->u
.param
->val
.nr
= sel
->v
.nr
;
744 if (sel
->u
.param
->nvalptr
)
746 *sel
->u
.param
->nvalptr
= sel
->u
.param
->val
.nr
;
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.
757 * If the value type is \ref POS_VALUE, the value of the child is simply
758 * copied to set the value of \p sel (the child subexpression should
759 * already have been evaluated by its root).
760 * If the value type is something else, the child is evaluated for the
761 * group \p g, and the value of the child is then copied.
762 * There should be only one child element.
764 * This function is used as gmx::SelectionTreeElement::evaluate for
765 * \ref SEL_SUBEXPRREF elements.
767 void _gmx_sel_evaluate_subexprref(gmx_sel_evaluate_t
* data
,
768 const gmx::SelectionTreeElementPointer
& sel
,
773 if (g
!= nullptr && sel
->child
->evaluate
!= nullptr)
775 sel
->child
->evaluate(data
, sel
->child
, g
);
777 const SelectionTreeElementPointer
& expr
= sel
->child
;
783 sel
->v
.nr
= expr
->v
.nr
;
784 memcpy(sel
->v
.u
.i
, expr
->v
.u
.i
, sel
->v
.nr
* sizeof(*sel
->v
.u
.i
));
788 sel
->v
.nr
= g
->isize
;
789 /* Extract the values corresponding to g */
790 for (i
= j
= 0; i
< g
->isize
; ++i
, ++j
)
792 while (sel
->child
->u
.cgrp
.index
[j
] < g
->index
[i
])
796 sel
->v
.u
.i
[i
] = expr
->v
.u
.i
[j
];
804 sel
->v
.nr
= expr
->v
.nr
;
805 memcpy(sel
->v
.u
.r
, expr
->v
.u
.r
, sel
->v
.nr
* sizeof(*sel
->v
.u
.r
));
809 sel
->v
.nr
= g
->isize
;
810 /* Extract the values corresponding to g */
811 for (i
= j
= 0; i
< g
->isize
; ++i
, ++j
)
813 while (sel
->child
->u
.cgrp
.index
[j
] < g
->index
[i
])
817 sel
->v
.u
.r
[i
] = expr
->v
.u
.r
[j
];
825 sel
->v
.nr
= expr
->v
.nr
;
826 memcpy(sel
->v
.u
.s
, expr
->v
.u
.s
, sel
->v
.nr
* sizeof(*sel
->v
.u
.s
));
830 sel
->v
.nr
= g
->isize
;
831 /* Extract the values corresponding to g */
832 for (i
= j
= 0; i
< g
->isize
; ++i
, ++j
)
834 while (sel
->child
->u
.cgrp
.index
[j
] < g
->index
[i
])
838 sel
->v
.u
.s
[i
] = expr
->v
.u
.s
[j
];
844 /* Currently, there is no need to do anything fancy here,
845 * but some future extensions may need a more flexible
847 gmx_ana_pos_copy(sel
->v
.u
.p
, expr
->v
.u
.p
, false);
853 gmx_ana_index_copy(sel
->v
.u
.g
, expr
->v
.u
.g
, false);
857 gmx_ana_index_intersection(sel
->v
.u
.g
, expr
->v
.u
.g
, g
);
861 default: /* should not be reached */
862 GMX_THROW(gmx::InternalError("Invalid subexpression reference type"));
864 /* Store the number of values if needed */
867 sel
->u
.param
->val
.nr
= sel
->v
.nr
;
868 if (sel
->u
.param
->nvalptr
)
870 *sel
->u
.param
->nvalptr
= sel
->u
.param
->val
.nr
;
875 /********************************************************************
876 * METHOD EXPRESSION EVALUATION
877 ********************************************************************/
880 * \param[in] data Data for the current frame.
881 * \param[in] sel Selection element being evaluated.
882 * \param[in] g Group for which \p sel should be evaluated.
883 * \returns 0 on success, a non-zero error code on error.
885 * Evaluates each child of a \ref SEL_EXPRESSION element.
886 * The value of \p sel is not touched.
888 * This function is not used as gmx::SelectionTreeElement::evaluate,
889 * but is used internally.
891 void _gmx_sel_evaluate_method_params(gmx_sel_evaluate_t
* data
,
892 const gmx::SelectionTreeElementPointer
& sel
,
895 SelectionTreeElementPointer child
= sel
->child
;
898 if (child
->evaluate
&& !(child
->flags
& SEL_EVALFRAME
))
900 if (child
->flags
& SEL_ATOMVAL
)
902 child
->evaluate(data
, child
, g
);
906 child
->flags
|= SEL_EVALFRAME
;
907 child
->evaluate(data
, child
, nullptr);
915 * \param[in] data Data for the current frame.
916 * \param[in] sel Selection element being evaluated.
917 * \param[in] g Group for which \p sel should be evaluated.
918 * \returns 0 on success, a non-zero error code on error.
920 * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
921 * to evaluate any parameter values.
922 * If this is the first time this expression is evaluated for
923 * the frame, sel_framefunc() callback is called if one is provided.
924 * If a reference position calculation has been initialized for this element,
925 * the positions are also updated, and sel_updatefunc_pos() is used to
926 * evaluate the value. Otherwise, sel_updatefunc() is used.
928 * This function is used as gmx::SelectionTreeElement::evaluate for
929 * \ref SEL_EXPRESSION elements.
931 void _gmx_sel_evaluate_method(gmx_sel_evaluate_t
* data
,
932 const gmx::SelectionTreeElementPointer
& sel
,
935 _gmx_sel_evaluate_method_params(data
, sel
, g
);
936 gmx::SelMethodEvalContext
context(data
->top
, data
->fr
, data
->pbc
);
937 if (sel
->flags
& SEL_INITFRAME
)
939 sel
->flags
&= ~SEL_INITFRAME
;
940 sel
->u
.expr
.method
->init_frame(context
, sel
->u
.expr
.mdata
);
944 gmx_ana_poscalc_update(sel
->u
.expr
.pc
, sel
->u
.expr
.pos
, g
, data
->fr
, data
->pbc
);
945 sel
->u
.expr
.method
->pupdate(context
, sel
->u
.expr
.pos
, &sel
->v
, sel
->u
.expr
.mdata
);
946 if ((sel
->flags
& SEL_ATOMVAL
) && sel
->v
.nr
< g
->isize
)
951 expandValueForPositions(sel
->v
.u
.r
, &sel
->v
.nr
, sel
->u
.expr
.pos
);
954 GMX_RELEASE_ASSERT(false,
955 "Unimplemented value type for position update method");
961 sel
->u
.expr
.method
->update(context
, g
, &sel
->v
, sel
->u
.expr
.mdata
);
966 * \param[in] data Data for the current frame.
967 * \param[in] sel Selection element being evaluated.
968 * \param[in] g Group for which \p sel should be evaluated.
969 * \returns 0 on success, a non-zero error code on error.
971 * Evaluates all child selections (using _gmx_sel_evaluate_method_params())
972 * to evaluate any parameter values.
973 * If this is the first time this expression is evaluated for
974 * the frame, sel_framefunc() callback is called if one is provided.
975 * The modifier is then evaluated using sel_updatefunc_pos().
977 * This function is used as gmx::SelectionTreeElement::evaluate for
978 * \ref SEL_MODIFIER elements.
980 void _gmx_sel_evaluate_modifier(gmx_sel_evaluate_t
* data
,
981 const gmx::SelectionTreeElementPointer
& sel
,
984 _gmx_sel_evaluate_method_params(data
, sel
, g
);
985 gmx::SelMethodEvalContext
context(data
->top
, data
->fr
, data
->pbc
);
986 if (sel
->flags
& SEL_INITFRAME
)
988 sel
->flags
&= ~SEL_INITFRAME
;
989 sel
->u
.expr
.method
->init_frame(context
, sel
->u
.expr
.mdata
);
991 if (sel
->child
&& sel
->child
->v
.type
!= POS_VALUE
)
993 GMX_THROW(gmx::NotImplementedError("Non-position valued modifiers not implemented"));
995 sel
->u
.expr
.method
->pupdate(context
, nullptr, &sel
->v
, sel
->u
.expr
.mdata
);
999 /********************************************************************
1000 * BOOLEAN EXPRESSION EVALUATION
1001 ********************************************************************/
1004 * \param[in] data Data for the current frame.
1005 * \param[in] sel Selection element being evaluated.
1006 * \param[in] g Group for which \p sel should be evaluated.
1007 * \returns 0 on success, a non-zero error code on error.
1009 * Evaluates the child element (there should be only one) in the group
1010 * \p g, and then sets the value of \p sel to the complement of the
1013 * This function is used as gmx::SelectionTreeElement::evaluate for
1014 * \ref SEL_BOOLEAN elements with \ref BOOL_NOT.
1016 void _gmx_sel_evaluate_not(gmx_sel_evaluate_t
* data
,
1017 const gmx::SelectionTreeElementPointer
& sel
,
1020 MempoolSelelemReserver
reserver(sel
->child
, g
->isize
);
1021 sel
->child
->evaluate(data
, sel
->child
, g
);
1022 gmx_ana_index_difference(sel
->v
.u
.g
, g
, sel
->child
->v
.u
.g
);
1026 * \param[in] data Data for the current frame.
1027 * \param[in] sel Selection element being evaluated.
1028 * \param[in] g Group for which \p sel should be evaluated.
1029 * \returns 0 on success, a non-zero error code on error.
1031 * Short-circuiting evaluation of logical AND expressions.
1033 * Starts by evaluating the first child element in the group \p g.
1034 * The each following child element is evaluated in the intersection
1035 * of all the previous values until all children have been evaluated
1036 * or the intersection becomes empty.
1037 * The value of \p sel is set to the intersection of all the (evaluated)
1040 * If the first child does not have an evaluation function, it is skipped
1041 * and the evaluation is started at the second child.
1042 * This happens if the first child is a constant expression and during
1043 * compilation it was detected that the evaluation group is always a subset
1044 * of the constant group
1045 * (currently, the compiler never detects this).
1047 * This function is used as gmx::SelectionTreeElement::evaluate for
1048 * \ref SEL_BOOLEAN elements with \ref BOOL_AND.
1050 void _gmx_sel_evaluate_and(gmx_sel_evaluate_t
* data
,
1051 const gmx::SelectionTreeElementPointer
& sel
,
1054 SelectionTreeElementPointer child
= sel
->child
;
1055 /* Skip the first child if it does not have an evaluation function. */
1056 if (!child
->evaluate
)
1058 child
= child
->next
;
1060 /* Evaluate the first child */
1062 MempoolSelelemReserver
reserver(child
, g
->isize
);
1063 child
->evaluate(data
, child
, g
);
1064 gmx_ana_index_copy(sel
->v
.u
.g
, child
->v
.u
.g
, false);
1066 child
= child
->next
;
1067 while (child
&& sel
->v
.u
.g
->isize
> 0)
1069 MempoolSelelemReserver
reserver(child
, sel
->v
.u
.g
->isize
);
1070 child
->evaluate(data
, child
, sel
->v
.u
.g
);
1071 gmx_ana_index_intersection(sel
->v
.u
.g
, sel
->v
.u
.g
, child
->v
.u
.g
);
1072 child
= child
->next
;
1077 * \param[in] data Data for the current frame.
1078 * \param[in] sel Selection element being evaluated.
1079 * \param[in] g Group for which \p sel should be evaluated.
1080 * \returns 0 on success, a non-zero error code on error.
1082 * Short-circuiting evaluation of logical OR expressions.
1084 * Starts by evaluating the first child element in the group \p g.
1085 * For each subsequent child, finds the part of \p g that is not
1086 * included the value of any previous child, and evaluates the child
1087 * in that group until the last child is evaluated or all of \p g
1088 * is included in some child value.
1089 * The value of \p sel is set to the union of all the (evaluated)
1092 * If the first child does not have an evaluation function, its value is
1093 * used without evaluation.
1094 * This happens if the first child is a constant expression, the selection
1095 * has been compiled, and the evaluation group is the same for each frame.
1096 * In this case, the compiler has taken care of that the child value is a
1097 * subset of \p g, making it unnecessary to evaluate it.
1099 * This function is used as gmx::SelectionTreeElement::evaluate for
1100 * \ref SEL_BOOLEAN elements with \ref BOOL_OR.
1102 void _gmx_sel_evaluate_or(gmx_sel_evaluate_t
* data
,
1103 const gmx::SelectionTreeElementPointer
& sel
,
1106 gmx_ana_index_t tmp
, tmp2
;
1108 SelectionTreeElementPointer child
= sel
->child
;
1109 if (child
->evaluate
)
1111 MempoolSelelemReserver
reserver(child
, g
->isize
);
1112 child
->evaluate(data
, child
, g
);
1113 gmx_ana_index_partition(sel
->v
.u
.g
, &tmp
, g
, child
->v
.u
.g
);
1117 gmx_ana_index_partition(sel
->v
.u
.g
, &tmp
, g
, child
->v
.u
.g
);
1119 child
= child
->next
;
1120 while (child
&& tmp
.isize
> 0)
1123 MempoolSelelemReserver
reserver(child
, tmp
.isize
);
1124 child
->evaluate(data
, child
, &tmp
);
1125 gmx_ana_index_partition(&tmp
, &tmp2
, &tmp
, child
->v
.u
.g
);
1127 sel
->v
.u
.g
->isize
+= tmp
.isize
;
1128 tmp
.isize
= tmp2
.isize
;
1129 tmp
.index
= tmp2
.index
;
1130 child
= child
->next
;
1132 gmx_ana_index_sort(sel
->v
.u
.g
);
1136 /********************************************************************
1137 * ARITHMETIC EVALUATION
1138 ********************************************************************/
1141 * \param[in] data Data for the current frame.
1142 * \param[in] sel Selection element being evaluated.
1143 * \param[in] g Group for which \p sel should be evaluated.
1144 * \returns 0 on success, a non-zero error code on error.
1146 void _gmx_sel_evaluate_arithmetic(gmx_sel_evaluate_t
* data
,
1147 const gmx::SelectionTreeElementPointer
& sel
,
1151 real lval
, rval
= 0., val
= 0.;
1153 const SelectionTreeElementPointer
& left
= sel
->child
;
1154 const SelectionTreeElementPointer
& right
= left
->next
;
1156 SelelemTemporaryValueAssigner assigner
;
1157 MempoolSelelemReserver reserver
;
1160 assigner
.assign(left
, *sel
);
1163 reserver
.reserve(right
, g
->isize
);
1166 else if (right
&& right
->mempool
)
1168 assigner
.assign(right
, *sel
);
1170 _gmx_sel_evaluate_children(data
, sel
, g
);
1172 n
= (sel
->flags
& SEL_SINGLEVAL
) ? 1 : g
->isize
;
1175 bool bArithNeg
= (sel
->u
.arith
.type
== ARITH_NEG
);
1176 GMX_ASSERT(right
|| bArithNeg
, "Right operand cannot be null except for negations");
1177 for (i
= i1
= i2
= 0; i
< n
; ++i
)
1179 lval
= left
->v
.u
.r
[i1
];
1182 rval
= right
->v
.u
.r
[i2
];
1184 switch (sel
->u
.arith
.type
)
1186 case ARITH_PLUS
: val
= lval
+ rval
; break;
1187 case ARITH_MINUS
: val
= lval
- rval
; break;
1188 case ARITH_NEG
: val
= -lval
; break;
1189 case ARITH_MULT
: val
= lval
* rval
; break;
1190 case ARITH_DIV
: val
= lval
/ rval
; break;
1191 case ARITH_EXP
: val
= pow(lval
, rval
); break;
1193 sel
->v
.u
.r
[i
] = val
;
1194 if (!(left
->flags
& SEL_SINGLEVAL
))
1198 if (!bArithNeg
&& !(right
->flags
& SEL_SINGLEVAL
))