Merge "Merge release-4-6 into master"
[alexxy/gromacs.git] / src / gromacs / selection / selection.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6  * others, as listed in the AUTHORS file in the top-level source
7  * directory and at http://www.gromacs.org.
8  *
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.
13  *
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.
18  *
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.
23  *
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.
31  *
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.
34  */
35 /*! \file
36  * \brief
37  * Declares gmx::Selection and supporting classes.
38  *
39  * \author Teemu Murtola <teemu.murtola@gmail.com>
40  * \inpublicapi
41  * \ingroup module_selection
42  */
43 #ifndef GMX_SELECTION_SELECTION_H
44 #define GMX_SELECTION_SELECTION_H
45
46 #include <string>
47 #include <vector>
48
49 #include "../legacyheaders/typedefs.h"
50
51 #include "../utility/arrayref.h"
52 #include "../utility/common.h"
53 #include "../utility/gmxassert.h"
54
55 #include "position.h"
56 #include "indexutil.h"
57 #include "selectionenums.h"
58
59 namespace gmx
60 {
61
62 class SelectionOptionStorage;
63 class SelectionTreeElement;
64
65 class Selection;
66 class SelectionPosition;
67
68 //! Container of selections used in public selection interfaces.
69 typedef std::vector<Selection> SelectionList;
70
71 namespace internal
72 {
73
74 /*! \internal \brief
75  * Internal data for a single selection.
76  *
77  * This class is internal to the selection module, but resides in a public
78  * header because of efficiency reasons: it allows frequently used access
79  * methods in \ref Selection to be inlined.
80  *
81  * Methods in this class do not throw unless otherwise specified.
82  *
83  * \ingroup module_selection
84  */
85 class SelectionData
86 {
87     public:
88         /*! \brief
89          * Creates a new selection object.
90          *
91          * \param[in] elem   Root of the evaluation tree for this selection.
92          * \param[in] selstr String that was parsed to produce this selection.
93          * \throws    std::bad_alloc if out of memory.
94          */
95         SelectionData(SelectionTreeElement *elem, const char *selstr);
96         ~SelectionData();
97
98         //! Returns the string that was parsed to produce this selection.
99         const char *selectionText() const { return selectionText_.c_str(); }
100         //! Returns true if the size of the selection (posCount()) is dynamic.
101         bool isDynamic() const { return bDynamic_; }
102         //! Number of positions in the selection.
103         int posCount() const { return rawPositions_.nr; }
104         //! Returns the root of the evaluation tree for this selection.
105         SelectionTreeElement &rootElement() { return rootElement_; }
106
107         //! Returns whether the covered fraction can change between frames.
108         bool isCoveredFractionDynamic() const { return bDynamicCoveredFraction_; }
109
110         //! Returns true if the given flag is set.
111         bool hasFlag(SelectionFlag flag) const { return flags_.test(flag); }
112         //! Sets the flags for this selection.
113         void setFlags(SelectionFlags flags) { flags_ = flags; }
114
115         //! \copydoc Selection::initCoveredFraction()
116         bool initCoveredFraction(e_coverfrac_t type);
117
118         /*! \brief
119          * Computes total masses and charges for all selection positions.
120          *
121          * \param[in] top   Topology information.
122          * \throws    std::bad_alloc if out of memory.
123          *
124          * For dynamic selections, the values need to be updated after each
125          * evaluation with refreshMassesAndCharges().
126          * This is done by SelectionEvaluator.
127          *
128          * This function is called by SelectionCompiler.
129          *
130          * Strong exception safety guarantee.
131          */
132         void initializeMassesAndCharges(const t_topology *top);
133         /*! \brief
134          * Updates masses and charges after dynamic selection has been
135          * evaluated.
136          *
137          * \param[in] top   Topology information.
138          *
139          * Called by SelectionEvaluator.
140          */
141         void refreshMassesAndCharges(const t_topology *top);
142         /*! \brief
143          * Updates the covered fraction after a selection has been evaluated.
144          *
145          * Called by SelectionEvaluator.
146          */
147         void updateCoveredFractionForFrame();
148         /*! \brief
149          * Computes average covered fraction after all frames have been evaluated.
150          *
151          * \param[in] nframes  Number of frames that have been evaluated.
152          *
153          * \p nframes should be equal to the number of calls to
154          * updateCoveredFractionForFrame().
155          * Called by SelectionEvaluator::evaluateFinal().
156          */
157         void computeAverageCoveredFraction(int nframes);
158         /*! \brief
159          * Restores position information to state it was in after compilation.
160          *
161          * \param[in] top   Topology information.
162          *
163          * Depends on SelectionCompiler storing the original atoms in the
164          * \a rootElement_ object.
165          * Called by SelectionEvaluator::evaluateFinal().
166          */
167         void restoreOriginalPositions(const t_topology *top);
168
169     private:
170         //! Name of the selection.
171         std::string               name_;
172         //! The actual selection string.
173         std::string               selectionText_;
174         //! Low-level representation of selected positions.
175         gmx_ana_pos_t             rawPositions_;
176         //! Total masses for the current positions.
177         std::vector<real>         posMass_;
178         //! Total charges for the current positions.
179         std::vector<real>         posCharge_;
180         SelectionFlags            flags_;
181         //! Root of the selection evaluation tree.
182         SelectionTreeElement     &rootElement_;
183         //! Type of the covered fraction.
184         e_coverfrac_t             coveredFractionType_;
185         //! Covered fraction of the selection for the current frame.
186         real                      coveredFraction_;
187         //! The average covered fraction (over the trajectory).
188         real                      averageCoveredFraction_;
189         //! true if the value can change as a function of time.
190         bool                      bDynamic_;
191         //! true if the covered fraction depends on the frame.
192         bool                      bDynamicCoveredFraction_;
193
194         /*! \brief
195          * Needed to wrap access to information.
196          */
197         friend class gmx::Selection;
198         /*! \brief
199          * Needed for proper access to position information.
200          */
201         friend class gmx::SelectionPosition;
202
203         GMX_DISALLOW_COPY_AND_ASSIGN(SelectionData);
204 };
205
206 }   // namespace internal
207
208 /*! \brief
209  * Provides access to a single selection.
210  *
211  * This class provides a public interface for accessing selection information.
212  * General information about the selection can be accessed with methods name(),
213  * selectionText(), isDynamic(), and type().  The first three can be accessed
214  * any time after the selection has been parsed, and type() can be accessed
215  * after the selection has been compiled.
216  *
217  * There are a few methods that can be used to change the behavior of the
218  * selection.  setEvaluateVelocities() and setEvaluateForces() can be called
219  * before the selection is compiled to request evaluation of velocities and/or
220  * forces in addition to coordinates.
221  *
222  * Each selection is made of a set of positions.  Each position has associated
223  * coordinates, and possibly velocities and forces if they have been requested
224  * and are available.  It also has a set of atoms associated with it; typically
225  * the coordinates are the center-of-mass or center-of-geometry coordinates for
226  * that set of atoms.  To access the number of positions in the selection, use
227  * posCount().  To access individual positions, use position().
228  * See SelectionPosition for details of how to use individual positions.
229  * setOriginalId() can be used to adjust the return value of
230  * SelectionPosition::mappedId(); see that method for details.
231  *
232  * It is also possible to access the list of atoms that make up all the
233  * positions directly: atomCount() returns the total number of atoms in the
234  * selection and atomIndices() an array of their indices.
235  * Similarly, it is possible to access the coordinates and other properties
236  * of the positions as continuous arrays through coordinates(), velocities(),
237  * forces(), masses(), charges(), refIds(), and mappedIds().
238  *
239  * Both positions and atoms can be accessed after the selection has been
240  * compiled.  For dynamic selections, the return values of these methods change
241  * after each evaluation to reflect the situation for the current frame.
242  * Before any frame has been evaluated, these methods return the maximal set
243  * to which the selection can evaluate.
244  *
245  * There are two possible modes for how positions for dynamic selections are
246  * handled.  In the default mode, posCount() can change, and for each frame,
247  * only the positions that are selected in that frame can be accessed.  In a
248  * masked mode, posCount() remains constant, i.e., the positions are always
249  * evaluated for the maximal set, and SelectionPosition::selected() is used to
250  * determine whether a position is selected for a frame.  The masked mode can
251  * be requested with SelectionOption::dynamicMask().
252  *
253  * The class also provides methods for printing out information: printInfo()
254  * and printDebugInfo().  These are mainly for internal use by Gromacs.
255  *
256  * This class works like a pointer type: copying and assignment is lightweight,
257  * and all copies work interchangeably, accessing the same internal data.
258  *
259  * Methods in this class do not throw.
260  *
261  * \see SelectionPosition
262  *
263  * \inpublicapi
264  * \ingroup module_selection
265  */
266 class Selection
267 {
268     public:
269         /*! \brief
270          * Creates a selection wrapper that has no associated selection.
271          *
272          * Any attempt to call methods in the object before a selection is
273          * assigned results in undefined behavior.
274          */
275         Selection() : sel_(NULL) {}
276         /*! \brief
277          * Creates a new selection object.
278          *
279          * \param  sel  Selection data to wrap.
280          *
281          * Only for internal use by the selection module.
282          */
283         explicit Selection(internal::SelectionData *sel) : sel_(sel) {}
284
285         //! Returns the name of the selection.
286         const char *name() const  { return data().name_.c_str(); }
287         //! Returns the string that was parsed to produce this selection.
288         const char *selectionText() const { return data().selectionText(); }
289         //! Returns true if the size of the selection (posCount()) is dynamic.
290         bool isDynamic() const { return data().isDynamic(); }
291         //! Returns the type of positions in the selection.
292         e_index_t type() const { return data().rawPositions_.m.type; }
293
294         //! Total number of atoms in the selection.
295         int atomCount() const
296         {
297             return data().rawPositions_.g != NULL ? data().rawPositions_.g->isize : 0;
298         }
299         //! Returns atom indices of all atoms in the selection.
300         ConstArrayRef<int> atomIndices() const
301         {
302             if (data().rawPositions_.g == NULL)
303             {
304                 return ConstArrayRef<int>();
305             }
306             return ConstArrayRef<int>(data().rawPositions_.g->index,
307                                       data().rawPositions_.g->isize);
308         }
309         //! Number of positions in the selection.
310         int posCount() const { return data().posCount(); }
311         //! Access a single position.
312         SelectionPosition position(int i) const;
313         //! Returns coordinates for this selection as a continuous array.
314         ConstArrayRef<rvec> coordinates() const
315         {
316             return ConstArrayRef<rvec>(data().rawPositions_.x, posCount());
317         }
318         //! Returns whether velocities are available for this selection.
319         bool hasVelocities() const { return data().rawPositions_.v != NULL; }
320         /*! \brief
321          * Returns velocities for this selection as a continuous array.
322          *
323          * Must not be called if hasVelocities() returns false.
324          */
325         ConstArrayRef<rvec> velocities() const
326         {
327             GMX_ASSERT(hasVelocities(), "Velocities accessed, but unavailable");
328             return ConstArrayRef<rvec>(data().rawPositions_.v, posCount());
329         }
330         //! Returns whether forces are available for this selection.
331         bool hasForces() const { return sel_->rawPositions_.f != NULL; }
332         /*! \brief
333          * Returns forces for this selection as a continuous array.
334          *
335          * Must not be called if hasForces() returns false.
336          */
337         ConstArrayRef<rvec> forces() const
338         {
339             GMX_ASSERT(hasForces(), "Forces accessed, but unavailable");
340             return ConstArrayRef<rvec>(data().rawPositions_.f, posCount());
341         }
342         //! Returns masses for this selection as a continuous array.
343         ConstArrayRef<real> masses() const
344         {
345             // posMass_ may have more entries than posCount() in the case of
346             // dynamic selections that don't have a topology
347             // (and thus the masses and charges are fixed).
348             GMX_ASSERT(data().posMass_.size() >= static_cast<size_t>(posCount()),
349                        "Internal inconsistency");
350             return ConstArrayRef<real>(data().posMass_.begin(),
351                                        data().posMass_.begin() + posCount());
352         }
353         //! Returns charges for this selection as a continuous array.
354         ConstArrayRef<real> charges() const
355         {
356             // posCharge_ may have more entries than posCount() in the case of
357             // dynamic selections that don't have a topology
358             // (and thus the masses and charges are fixed).
359             GMX_ASSERT(data().posCharge_.size() >= static_cast<size_t>(posCount()),
360                        "Internal inconsistency");
361             return ConstArrayRef<real>(data().posCharge_.begin(),
362                                        data().posCharge_.begin() + posCount());
363         }
364         /*! \brief
365          * Returns reference IDs for this selection as a continuous array.
366          *
367          * \see SelectionPosition::refId()
368          */
369         ConstArrayRef<int> refIds() const
370         {
371             return ConstArrayRef<int>(data().rawPositions_.m.refid, posCount());
372         }
373         /*! \brief
374          * Returns mapped IDs for this selection as a continuous array.
375          *
376          * \see SelectionPosition::mappedId()
377          */
378         ConstArrayRef<int> mappedIds() const
379         {
380             return ConstArrayRef<int>(data().rawPositions_.m.mapid, posCount());
381         }
382
383         //! Deprecated method for direct access to position data.
384         const gmx_ana_pos_t *positions() const { return &data().rawPositions_; }
385
386         //! Returns whether the covered fraction can change between frames.
387         bool isCoveredFractionDynamic() const { return data().isCoveredFractionDynamic(); }
388         //! Returns the covered fraction for the current frame.
389         real coveredFraction() const { return data().coveredFraction_; }
390
391         /*! \brief
392          * Initializes information about covered fractions.
393          *
394          * \param[in] type Type of covered fraction required.
395          * \returns   true if the covered fraction can be calculated for the
396          *      selection.
397          */
398         bool initCoveredFraction(e_coverfrac_t type)
399         {
400             return data().initCoveredFraction(type);
401         }
402         /*! \brief
403          * Sets whether this selection evaluates velocities for positions.
404          *
405          * \param[in] bEnabled  If true, velocities are evaluated.
406          *
407          * If you request the evaluation, but then evaluate the selection for
408          * a frame that does not contain velocity information, results are
409          * undefined.
410          *
411          * \todo
412          * Implement it such that in the above case, hasVelocities() will
413          * return false for such frames.
414          *
415          * Does not throw.
416          */
417         void setEvaluateVelocities(bool bEnabled)
418         {
419             data().flags_.set(efSelection_EvaluateVelocities, bEnabled);
420         }
421         /*! \brief
422          * Sets whether this selection evaluates forces for positions.
423          *
424          * \param[in] bEnabled  If true, forces are evaluated.
425          *
426          * If you request the evaluation, but then evaluate the selection for
427          * a frame that does not contain force information, results are
428          * undefined.
429          *
430          * Does not throw.
431          */
432         void setEvaluateForces(bool bEnabled)
433         {
434             data().flags_.set(efSelection_EvaluateForces, bEnabled);
435         }
436
437         /*! \brief
438          * Sets the ID for the \p i'th position for use with
439          * SelectionPosition::mappedId().
440          *
441          * \param[in] i  Zero-based index
442          * \param[in] id Identifier to set.
443          *
444          * This method is not part of SelectionPosition because that interface
445          * only provides access to const data by design.
446          *
447          * This method can only be called after compilation, before the
448          * selection has been evaluated for any frame.
449          *
450          * \see SelectionPosition::mappedId()
451          */
452         void setOriginalId(int i, int id) { data().rawPositions_.m.orgid[i] = id; }
453
454         /*! \brief
455          * Prints out one-line description of the selection.
456          *
457          * \param[in] fp      Where to print the information.
458          *
459          * The output contains the name of the selection, the number of atoms
460          * and the number of positions, and indication of whether the selection
461          * is dynamic.
462          */
463         void printInfo(FILE *fp) const;
464         /*! \brief
465          * Prints out extended information about the selection for debugging.
466          *
467          * \param[in] fp      Where to print the information.
468          * \param[in] nmaxind Maximum number of values to print in lists
469          *      (-1 = print all).
470          */
471         void printDebugInfo(FILE *fp, int nmaxind) const;
472
473     private:
474         internal::SelectionData &data()
475         {
476             GMX_ASSERT(sel_ != NULL,
477                        "Attempted to access uninitialized selection");
478             return *sel_;
479         }
480         const internal::SelectionData &data() const
481         {
482             GMX_ASSERT(sel_ != NULL,
483                        "Attempted to access uninitialized selection");
484             return *sel_;
485         }
486
487         /*! \brief
488          * Pointer to internal data for the selection.
489          *
490          * The memory for this object is managed by a SelectionCollection
491          * object, and the \ref Selection class simply provides a public
492          * interface for accessing the data.
493          */
494         internal::SelectionData *sel_;
495
496         /*! \brief
497          * Needed to access the data to adjust flags.
498          */
499         friend class SelectionOptionStorage;
500 };
501
502 /*! \brief
503  * Provides access to information about a single selected position.
504  *
505  * Each position has associated coordinates, and possibly velocities and forces
506  * if they have been requested and are available.  It also has a set of atoms
507  * associated with it; typically the coordinates are the center-of-mass or
508  * center-of-geometry coordinates for that set of atoms.  It is possible that
509  * there are not atoms associated if the selection has been provided as a fixed
510  * position.
511  *
512  * After the selection has been compiled, but not yet evaluated, the contents
513  * of the coordinate, velocity and force vectors are undefined.
514  *
515  * Default copy constructor and assignment operators are used, and work as
516  * intended: the copy references the same position and works identically.
517  *
518  * Methods in this class do not throw.
519  *
520  * \see Selection
521  *
522  * \inpublicapi
523  * \ingroup module_selection
524  */
525 class SelectionPosition
526 {
527     public:
528         /*! \brief
529          * Constructs a wrapper object for given selection position.
530          *
531          * \param[in] sel    Selection from which the position is wrapped.
532          * \param[in] index  Zero-based index of the position to wrap.
533          *
534          * Asserts if \p index is out of range.
535          *
536          * Only for internal use of the library.  To obtain a SelectionPosition
537          * object in other code, use Selection::position().
538          */
539         SelectionPosition(const internal::SelectionData &sel, int index)
540             : sel_(&sel), i_(index)
541         {
542             GMX_ASSERT(index >= 0 && index < sel.posCount(),
543                        "Invalid selection position index");
544         }
545
546         /*! \brief
547          * Returns type of this position.
548          *
549          * Currently always returns the same as Selection::type().
550          */
551         e_index_t type() const { return sel_->rawPositions_.m.type; }
552         //! Returns coordinates for this position.
553         const rvec &x() const
554         {
555             return sel_->rawPositions_.x[i_];
556         }
557         /*! \brief
558          * Returns velocity for this position.
559          *
560          * Must not be called if Selection::hasVelocities() returns false.
561          */
562         const rvec &v() const
563         {
564             GMX_ASSERT(sel_->rawPositions_.v != NULL,
565                        "Velocities accessed, but unavailable");
566             return sel_->rawPositions_.v[i_];
567         }
568         /*! \brief
569          * Returns force for this position.
570          *
571          * Must not be called if Selection::hasForces() returns false.
572          */
573         const rvec &f() const
574         {
575             GMX_ASSERT(sel_->rawPositions_.f != NULL,
576                        "Velocities accessed, but unavailable");
577             return sel_->rawPositions_.f[i_];
578         }
579         /*! \brief
580          * Returns total mass for this position.
581          *
582          * Returns the total mass of atoms that make up this position.
583          * If there are no atoms associated or masses are not available,
584          * returns unity.
585          */
586         real mass() const
587         {
588             return sel_->posMass_[i_];
589         }
590         /*! \brief
591          * Returns total charge for this position.
592          *
593          * Returns the sum of charges of atoms that make up this position.
594          * If there are no atoms associated or charges are not available,
595          * returns zero.
596          */
597         real charge() const
598         {
599             return sel_->posCharge_[i_];
600         }
601         //! Returns the number of atoms that make up this position.
602         int atomCount() const
603         {
604             return sel_->rawPositions_.m.mapb.index[i_ + 1]
605                    - sel_->rawPositions_.m.mapb.index[i_];
606         }
607         //! Return atom indices that make up this position.
608         ConstArrayRef<int> atomIndices() const
609         {
610             if (sel_->rawPositions_.g == NULL)
611             {
612                 return ConstArrayRef<int>();
613             }
614             int first = sel_->rawPositions_.m.mapb.index[i_];
615             return ConstArrayRef<int>(&sel_->rawPositions_.g->index[first],
616                                       atomCount());
617         }
618         /*! \brief
619          * Returns whether this position is selected in the current frame.
620          *
621          * The return value is equivalent to \c refid() == -1.  Returns always
622          * true if SelectionOption::dynamicMask() has not been set.
623          *
624          * \see refId()
625          */
626         bool selected() const
627         {
628             return refId() >= 0;
629         }
630         /*! \brief
631          * Returns reference ID for this position.
632          *
633          * For dynamic selections, this provides means to associate positions
634          * across frames.  After compilation, these IDs are consequently
635          * numbered starting from zero.  For each frame, the ID then reflects
636          * the location of the position in the original array of positions.
637          * If SelectionOption::dynamicMask() has been set for the parent
638          * selection, the IDs for positions not present in the current
639          * selection are set to -1, otherwise they are removed completely.
640          *
641          * Example:
642          * If a dynamic selection consists of at most three positions, after
643          * compilation refId() will return 0, 1, 2 for them, respectively.
644          * If for a particular frame, only the first and the third are present,
645          * refId() will return 0, 2.
646          * If SelectionOption::dynamicMask() has been set, all three positions
647          * can be accessed also for that frame and refId() will return 0, -1,
648          * 2.
649          */
650         int refId() const
651         {
652             return sel_->rawPositions_.m.refid[i_];
653         }
654         /*! \brief
655          * Returns mapped ID for this position.
656          *
657          * Returns ID of the position that corresponds to that set with
658          * Selection::setOriginalId().
659          *
660          * If for an array \c id, \c setOriginalId(i, id[i]) has been called
661          * for each \c i, then it always holds that
662          * \c mappedId()==id[refId()].
663          *
664          * Selection::setOriginalId() has not been called, the default values
665          * are dependent on type():
666          *  - ::INDEX_ATOM: atom indices
667          *  - ::INDEX_RES:  residue numbers
668          *  - ::INDEX_MOL:  molecule numbers
669          *  .
670          * All the default values are zero-based
671          */
672         int mappedId() const
673         {
674             return sel_->rawPositions_.m.mapid[i_];
675         }
676
677     private:
678         const internal::SelectionData  *sel_;
679         int                             i_;
680 };
681
682
683 inline SelectionPosition
684 Selection::position(int i) const
685 {
686     return SelectionPosition(data(), i);
687 }
688
689 } // namespace gmx
690
691 #endif