288bcea5e0d9ae5b36ade0f859c3a1d47b6c2445
[alexxy/gromacs.git] / src / gromacs / selection / selectioncollection.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2010,2011,2012,2013,2014, 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.
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 /*! \internal \file
36  * \brief
37  * Implements gmx::SelectionCollection.
38  *
39  * \author Teemu Murtola <teemu.murtola@gmail.com>
40  * \ingroup module_selection
41  */
42 #include "gmxpre.h"
43
44 #include "selectioncollection.h"
45
46 #include <cctype>
47 #include <cstdio>
48
49 #include <string>
50 #include <vector>
51
52 #include <boost/shared_ptr.hpp>
53
54 #include "gromacs/legacyheaders/oenv.h"
55
56 #include "gromacs/fileio/trx.h"
57 #include "gromacs/onlinehelp/helpmanager.h"
58 #include "gromacs/onlinehelp/helpwritercontext.h"
59 #include "gromacs/options/basicoptions.h"
60 #include "gromacs/options/options.h"
61 #include "gromacs/selection/selection.h"
62 #include "gromacs/topology/topology.h"
63 #include "gromacs/utility/exceptions.h"
64 #include "gromacs/utility/file.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/messagestringcollector.h"
67 #include "gromacs/utility/smalloc.h"
68 #include "gromacs/utility/stringutil.h"
69
70 #include "compiler.h"
71 #include "mempool.h"
72 #include "parser.h"
73 #include "poscalc.h"
74 #include "scanner.h"
75 #include "selection.h"
76 #include "selectioncollection-impl.h"
77 #include "selelem.h"
78 #include "selhelp.h"
79 #include "selmethod.h"
80 #include "symrec.h"
81
82 namespace gmx
83 {
84
85 /********************************************************************
86  * SelectionCollection::Impl
87  */
88
89 SelectionCollection::Impl::Impl()
90     : maxAtomIndex_(0), debugLevel_(0), bExternalGroupsSet_(false), grps_(NULL)
91 {
92     sc_.nvars     = 0;
93     sc_.varstrs   = NULL;
94     sc_.top       = NULL;
95     gmx_ana_index_clear(&sc_.gall);
96     sc_.mempool   = NULL;
97     sc_.symtab.reset(new SelectionParserSymbolTable);
98     gmx_ana_selmethod_register_defaults(sc_.symtab.get());
99 }
100
101
102 SelectionCollection::Impl::~Impl()
103 {
104     clearSymbolTable();
105     // The tree must be freed before the SelectionData objects, since the
106     // tree may hold references to the position data in SelectionData.
107     sc_.root.reset();
108     sc_.sel.clear();
109     for (int i = 0; i < sc_.nvars; ++i)
110     {
111         sfree(sc_.varstrs[i]);
112     }
113     sfree(sc_.varstrs);
114     gmx_ana_index_deinit(&sc_.gall);
115     if (sc_.mempool)
116     {
117         _gmx_sel_mempool_destroy(sc_.mempool);
118     }
119 }
120
121
122 void
123 SelectionCollection::Impl::clearSymbolTable()
124 {
125     sc_.symtab.reset();
126 }
127
128
129 namespace
130 {
131
132 /*! \brief
133  * Reads a single selection line from stdin.
134  *
135  * \param[in]  infile        File to read from (typically File::standardInput()).
136  * \param[in]  bInteractive  Whether to print interactive prompts.
137  * \param[out] line          The read line in stored here.
138  * \returns true if something was read, false if at end of input.
139  *
140  * Handles line continuation, reading also the continuing line(s) in one call.
141  */
142 bool promptLine(File *infile, bool bInteractive, std::string *line)
143 {
144     if (bInteractive)
145     {
146         fprintf(stderr, "> ");
147     }
148     if (!infile->readLineWithTrailingSpace(line))
149     {
150         return false;
151     }
152     while (endsWith(*line, "\\\n"))
153     {
154         line->resize(line->length() - 2);
155         if (bInteractive)
156         {
157             fprintf(stderr, "... ");
158         }
159         std::string buffer;
160         // Return value ignored, buffer remains empty and works correctly
161         // if there is nothing to read.
162         infile->readLineWithTrailingSpace(&buffer);
163         line->append(buffer);
164     }
165     if (endsWith(*line, "\n"))
166     {
167         line->resize(line->length() - 1);
168     }
169     else if (bInteractive)
170     {
171         fprintf(stderr, "\n");
172     }
173     return true;
174 }
175
176 /*! \brief
177  * Helper function for tokenizing the input and pushing them to the parser.
178  *
179  * \param     scanner       Tokenizer data structure.
180  * \param     parserState   Parser data structure.
181  * \param[in] bInteractive  Whether to operate in interactive mode.
182  *
183  * Repeatedly reads tokens using \p scanner and pushes them to the parser with
184  * \p parserState until there is no more input, or until enough input is given
185  * (only in interactive mode).
186  */
187 int runParserLoop(yyscan_t scanner, _gmx_sel_yypstate *parserState,
188                   bool bInteractive)
189 {
190     int status    = YYPUSH_MORE;
191     int prevToken = 0;
192     do
193     {
194         YYSTYPE value;
195         int     token = _gmx_sel_yylex(&value, scanner);
196         if (bInteractive)
197         {
198             if (token == 0)
199             {
200                 break;
201             }
202             // Empty commands cause the interactive parser to print out
203             // status information. This avoids producing those unnecessarily,
204             // e.g., from "resname RA;;".
205             if (prevToken == CMD_SEP && token == CMD_SEP)
206             {
207                 continue;
208             }
209             prevToken = token;
210         }
211         status = _gmx_sel_yypush_parse(parserState, token, &value, scanner);
212     }
213     while (status == YYPUSH_MORE);
214     _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
215     return status;
216 }
217
218 /*! \brief
219  * Print current status in response to empty line in interactive input.
220  *
221  * \param[in] sc             Selection collection data structure.
222  * \param[in] grps           Available index groups.
223  * \param[in] firstSelection Index of first selection from this interactive
224  *     session.
225  * \param[in] maxCount       Maximum number of selections.
226  * \param[in] context        Context to print for what the selections are for.
227  * \param[in] bFirst         Whether this is the header that is printed before
228  *     any user input.
229  *
230  * Prints the available index groups and currently provided selections.
231  */
232 void printCurrentStatus(gmx_ana_selcollection_t *sc, gmx_ana_indexgrps_t *grps,
233                         size_t firstSelection, int maxCount,
234                         const std::string &context, bool bFirst)
235 {
236     if (grps != NULL)
237     {
238         std::fprintf(stderr, "Available static index groups:\n");
239         gmx_ana_indexgrps_print(stderr, grps, 0);
240     }
241     std::fprintf(stderr, "Specify ");
242     if (maxCount < 0)
243     {
244         std::fprintf(stderr, "any number of selections");
245     }
246     else if (maxCount == 1)
247     {
248         std::fprintf(stderr, "a selection");
249     }
250     else
251     {
252         std::fprintf(stderr, "%d selections", maxCount);
253     }
254     std::fprintf(stderr, "%s%s:\n",
255                  context.empty() ? "" : " ", context.c_str());
256     std::fprintf(stderr,
257                  "(one per line, <enter> for status/groups, 'help' for help%s)\n",
258                  maxCount < 0 ? ", Ctrl-D to end" : "");
259     if (!bFirst && (sc->nvars > 0 || sc->sel.size() > firstSelection))
260     {
261         std::fprintf(stderr, "Currently provided selections:\n");
262         for (int i = 0; i < sc->nvars; ++i)
263         {
264             std::fprintf(stderr, "     %s\n", sc->varstrs[i]);
265         }
266         for (size_t i = firstSelection; i < sc->sel.size(); ++i)
267         {
268             std::fprintf(stderr, " %2d. %s\n",
269                          static_cast<int>(i - firstSelection + 1),
270                          sc->sel[i]->selectionText());
271         }
272         if (maxCount > 0)
273         {
274             const int remaining
275                 = maxCount - static_cast<int>(sc->sel.size() - firstSelection);
276             std::fprintf(stderr, "(%d more selection%s required)\n",
277                          remaining, remaining > 1 ? "s" : "");
278         }
279     }
280 }
281
282 /*! \brief
283  * Prints selection help in interactive selection input.
284  *
285  * \param[in] sc    Selection collection data structure.
286  * \param[in] line  Line of user input requesting help (starting with `help`).
287  *
288  * Initializes the selection help if not yet initialized, and finds the help
289  * topic based on words on the input line.
290  */
291 void printHelp(gmx_ana_selcollection_t *sc, const std::string &line)
292 {
293     if (sc->rootHelp.get() == NULL)
294     {
295         sc->rootHelp = createSelectionHelpTopic();
296     }
297     HelpWriterContext context(&File::standardError(),
298                               eHelpOutputFormat_Console);
299     HelpManager       manager(*sc->rootHelp, context);
300     try
301     {
302         std::vector<std::string>                 topic = splitString(line);
303         std::vector<std::string>::const_iterator value;
304         // First item in the list is the 'help' token.
305         for (value = topic.begin() + 1; value != topic.end(); ++value)
306         {
307             manager.enterTopic(*value);
308         }
309     }
310     catch (const InvalidInputError &ex)
311     {
312         fprintf(stderr, "%s\n", ex.what());
313         return;
314     }
315     manager.writeCurrentTopic();
316 }
317
318 /*! \brief
319  * Helper function that runs the parser once the tokenizer has been
320  * initialized.
321  *
322  * \param[in,out] scanner Scanner data structure.
323  * \param[in]     bStdIn  Whether to use a line-based reading
324  *      algorithm designed for interactive input.
325  * \param[in]     maxnr   Maximum number of selections to parse
326  *      (if -1, parse as many as provided by the user).
327  * \param[in]     context Context to print for what the selections are for.
328  * \returns       Vector of parsed selections.
329  * \throws        std::bad_alloc if out of memory.
330  * \throws        InvalidInputError if there is a parsing error.
331  *
332  * Used internally to implement parseFromStdin(), parseFromFile() and
333  * parseFromString().
334  */
335 SelectionList runParser(yyscan_t scanner, bool bStdIn, int maxnr,
336                         const std::string &context)
337 {
338     boost::shared_ptr<void>  scannerGuard(scanner, &_gmx_sel_free_lexer);
339     gmx_ana_selcollection_t *sc   = _gmx_sel_lexer_selcollection(scanner);
340     gmx_ana_indexgrps_t     *grps = _gmx_sel_lexer_indexgrps(scanner);
341
342     MessageStringCollector   errors;
343     _gmx_sel_set_lexer_error_reporter(scanner, &errors);
344
345     size_t oldCount = sc->sel.size();
346     bool   bOk      = false;
347     {
348         boost::shared_ptr<_gmx_sel_yypstate> parserState(
349                 _gmx_sel_yypstate_new(), &_gmx_sel_yypstate_delete);
350         if (bStdIn)
351         {
352             File       &stdinFile(File::standardInput());
353             const bool  bInteractive = _gmx_sel_is_lexer_interactive(scanner);
354             if (bInteractive)
355             {
356                 printCurrentStatus(sc, grps, oldCount, maxnr, context, true);
357             }
358             std::string line;
359             int         status;
360             while (promptLine(&stdinFile, bInteractive, &line))
361             {
362                 if (bInteractive)
363                 {
364                     line = stripString(line);
365                     if (line.empty())
366                     {
367                         printCurrentStatus(sc, grps, oldCount, maxnr, context, false);
368                         continue;
369                     }
370                     if (startsWith(line, "help")
371                         && (line[4] == 0 || std::isspace(line[4])))
372                     {
373                         printHelp(sc, line);
374                         continue;
375                     }
376                 }
377                 line.append("\n");
378                 _gmx_sel_set_lex_input_str(scanner, line.c_str());
379                 status = runParserLoop(scanner, parserState.get(), true);
380                 if (status != YYPUSH_MORE)
381                 {
382                     // TODO: Check if there is more input, and issue an
383                     // error/warning if some input was ignored.
384                     goto early_termination;
385                 }
386                 if (!errors.isEmpty() && bInteractive)
387                 {
388                     fprintf(stderr, "%s", errors.toString().c_str());
389                     errors.clear();
390                 }
391             }
392             status = _gmx_sel_yypush_parse(parserState.get(), 0, NULL,
393                                            scanner);
394             _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
395 early_termination:
396             bOk = (status == 0);
397         }
398         else
399         {
400             int status = runParserLoop(scanner, parserState.get(), false);
401             bOk = (status == 0);
402         }
403     }
404     scannerGuard.reset();
405     int nr = sc->sel.size() - oldCount;
406     if (maxnr > 0 && nr != maxnr)
407     {
408         bOk = false;
409         errors.append("Too few selections provided");
410     }
411
412     // TODO: Remove added selections from the collection if parsing failed?
413     if (!bOk || !errors.isEmpty())
414     {
415         GMX_ASSERT(!bOk && !errors.isEmpty(), "Inconsistent error reporting");
416         GMX_THROW(InvalidInputError(errors.toString()));
417     }
418
419     SelectionList                     result;
420     SelectionDataList::const_iterator i;
421     result.reserve(nr);
422     for (i = sc->sel.begin() + oldCount; i != sc->sel.end(); ++i)
423     {
424         result.push_back(Selection(i->get()));
425     }
426     return result;
427 }
428
429 /*! \brief
430  * Checks that index groups have valid atom indices.
431  *
432  * \param[in]    root    Root of selection tree to process.
433  * \param[in]    natoms  Maximum number of atoms that the selections are set
434  *     to evaluate.
435  * \param        errors  Object for reporting any error messages.
436  * \throws std::bad_alloc if out of memory.
437  *
438  * Recursively checks the selection tree for index groups.
439  * Each found group is checked that it only contains atom indices that match
440  * the topology/maximum number of atoms set for the selection collection.
441  * Any issues are reported to \p errors.
442  */
443 void checkExternalGroups(const SelectionTreeElementPointer &root,
444                          int                                natoms,
445                          ExceptionInitializer              *errors)
446 {
447     if (root->type == SEL_CONST && root->v.type == GROUP_VALUE)
448     {
449         try
450         {
451             root->checkIndexGroup(natoms);
452         }
453         catch (const UserInputError &)
454         {
455             errors->addCurrentExceptionAsNested();
456         }
457     }
458
459     SelectionTreeElementPointer child = root->child;
460     while (child)
461     {
462         checkExternalGroups(child, natoms, errors);
463         child = child->next;
464     }
465 }
466
467 }   // namespace
468
469
470 void SelectionCollection::Impl::resolveExternalGroups(
471         const SelectionTreeElementPointer &root,
472         ExceptionInitializer              *errors)
473 {
474
475     if (root->type == SEL_GROUPREF)
476     {
477         try
478         {
479             root->resolveIndexGroupReference(grps_, sc_.gall.isize);
480         }
481         catch (const UserInputError &)
482         {
483             errors->addCurrentExceptionAsNested();
484         }
485     }
486
487     SelectionTreeElementPointer child = root->child;
488     while (child)
489     {
490         resolveExternalGroups(child, errors);
491         root->flags |= (child->flags & SEL_UNSORTED);
492         child        = child->next;
493     }
494 }
495
496
497 /********************************************************************
498  * SelectionCollection
499  */
500
501 SelectionCollection::SelectionCollection()
502     : impl_(new Impl)
503 {
504 }
505
506
507 SelectionCollection::~SelectionCollection()
508 {
509 }
510
511
512 void
513 SelectionCollection::initOptions(Options *options)
514 {
515     const char * const debug_levels[]
516         = { "no", "basic", "compile", "eval", "full" };
517
518     bool bAllowNonAtomOutput = false;
519     SelectionDataList::const_iterator iter;
520     for (iter = impl_->sc_.sel.begin(); iter != impl_->sc_.sel.end(); ++iter)
521     {
522         const internal::SelectionData &sel = **iter;
523         if (!sel.hasFlag(efSelection_OnlyAtoms))
524         {
525             bAllowNonAtomOutput = true;
526         }
527     }
528
529     const char *const *postypes = PositionCalculationCollection::typeEnumValues;
530     options->addOption(StringOption("selrpos")
531                            .enumValueFromNullTerminatedArray(postypes)
532                            .store(&impl_->rpost_).defaultValue(postypes[0])
533                            .description("Selection reference positions"));
534     if (bAllowNonAtomOutput)
535     {
536         options->addOption(StringOption("seltype")
537                                .enumValueFromNullTerminatedArray(postypes)
538                                .store(&impl_->spost_).defaultValue(postypes[0])
539                                .description("Default selection output positions"));
540     }
541     else
542     {
543         impl_->spost_ = postypes[0];
544     }
545     GMX_RELEASE_ASSERT(impl_->debugLevel_ >= 0 && impl_->debugLevel_ <= 4,
546                        "Debug level out of range");
547     options->addOption(StringOption("seldebug").hidden(impl_->debugLevel_ == 0)
548                            .enumValue(debug_levels)
549                            .defaultValue(debug_levels[impl_->debugLevel_])
550                            .storeEnumIndex(&impl_->debugLevel_)
551                            .description("Print out selection trees for debugging"));
552 }
553
554
555 void
556 SelectionCollection::setReferencePosType(const char *type)
557 {
558     GMX_RELEASE_ASSERT(type != NULL, "Cannot assign NULL position type");
559     // Check that the type is valid, throw if it is not.
560     e_poscalc_t  dummytype;
561     int          dummyflags;
562     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
563     impl_->rpost_ = type;
564 }
565
566
567 void
568 SelectionCollection::setOutputPosType(const char *type)
569 {
570     GMX_RELEASE_ASSERT(type != NULL, "Cannot assign NULL position type");
571     // Check that the type is valid, throw if it is not.
572     e_poscalc_t  dummytype;
573     int          dummyflags;
574     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
575     impl_->spost_ = type;
576 }
577
578
579 void
580 SelectionCollection::setDebugLevel(int debugLevel)
581 {
582     impl_->debugLevel_ = debugLevel;
583 }
584
585
586 void
587 SelectionCollection::setTopology(t_topology *top, int natoms)
588 {
589     GMX_RELEASE_ASSERT(natoms > 0 || top != NULL,
590                        "The number of atoms must be given if there is no topology");
591     // Get the number of atoms from the topology if it is not given.
592     if (natoms <= 0)
593     {
594         natoms = top->atoms.nr;
595     }
596     if (impl_->bExternalGroupsSet_)
597     {
598         ExceptionInitializer        errors("Invalid index group references encountered");
599         SelectionTreeElementPointer root = impl_->sc_.root;
600         while (root)
601         {
602             checkExternalGroups(root, natoms, &errors);
603             root = root->next;
604         }
605         if (errors.hasNestedExceptions())
606         {
607             GMX_THROW(InconsistentInputError(errors));
608         }
609     }
610     gmx_ana_selcollection_t *sc = &impl_->sc_;
611     // Do this first, as it allocates memory, while the others don't throw.
612     gmx_ana_index_init_simple(&sc->gall, natoms);
613     sc->pcc.setTopology(top);
614     sc->top = top;
615 }
616
617
618 void
619 SelectionCollection::setIndexGroups(gmx_ana_indexgrps_t *grps)
620 {
621     GMX_RELEASE_ASSERT(grps == NULL || !impl_->bExternalGroupsSet_,
622                        "Can only set external groups once or clear them afterwards");
623     impl_->grps_               = grps;
624     impl_->bExternalGroupsSet_ = true;
625
626     ExceptionInitializer        errors("Invalid index group reference(s)");
627     SelectionTreeElementPointer root = impl_->sc_.root;
628     while (root)
629     {
630         impl_->resolveExternalGroups(root, &errors);
631         root->checkUnsortedAtoms(true, &errors);
632         root = root->next;
633     }
634     if (errors.hasNestedExceptions())
635     {
636         GMX_THROW(InconsistentInputError(errors));
637     }
638     for (size_t i = 0; i < impl_->sc_.sel.size(); ++i)
639     {
640         impl_->sc_.sel[i]->refreshName();
641     }
642 }
643
644
645 bool
646 SelectionCollection::requiresTopology() const
647 {
648     e_poscalc_t  type;
649     int          flags;
650
651     if (!impl_->rpost_.empty())
652     {
653         flags = 0;
654         // Should not throw, because has been checked earlier.
655         PositionCalculationCollection::typeFromEnum(impl_->rpost_.c_str(),
656                                                     &type, &flags);
657         if (type != POS_ATOM)
658         {
659             return true;
660         }
661     }
662     if (!impl_->spost_.empty())
663     {
664         flags = 0;
665         // Should not throw, because has been checked earlier.
666         PositionCalculationCollection::typeFromEnum(impl_->spost_.c_str(),
667                                                     &type, &flags);
668         if (type != POS_ATOM)
669         {
670             return true;
671         }
672     }
673
674     SelectionTreeElementPointer sel = impl_->sc_.root;
675     while (sel)
676     {
677         if (_gmx_selelem_requires_top(*sel))
678         {
679             return true;
680         }
681         sel = sel->next;
682     }
683     return false;
684 }
685
686
687 SelectionList
688 SelectionCollection::parseFromStdin(int nr, bool bInteractive,
689                                     const std::string &context)
690 {
691     yyscan_t scanner;
692
693     _gmx_sel_init_lexer(&scanner, &impl_->sc_, bInteractive, nr,
694                         impl_->bExternalGroupsSet_,
695                         impl_->grps_);
696     return runParser(scanner, true, nr, context);
697 }
698
699
700 SelectionList
701 SelectionCollection::parseFromFile(const std::string &filename)
702 {
703
704     try
705     {
706         yyscan_t scanner;
707         File     file(filename, "r");
708         // TODO: Exception-safe way of using the lexer.
709         _gmx_sel_init_lexer(&scanner, &impl_->sc_, false, -1,
710                             impl_->bExternalGroupsSet_,
711                             impl_->grps_);
712         _gmx_sel_set_lex_input_file(scanner, file.handle());
713         return runParser(scanner, false, -1, std::string());
714     }
715     catch (GromacsException &ex)
716     {
717         ex.prependContext(formatString(
718                                   "Error in parsing selections from file '%s'",
719                                   filename.c_str()));
720         throw;
721     }
722 }
723
724
725 SelectionList
726 SelectionCollection::parseFromString(const std::string &str)
727 {
728     yyscan_t scanner;
729
730     _gmx_sel_init_lexer(&scanner, &impl_->sc_, false, -1,
731                         impl_->bExternalGroupsSet_,
732                         impl_->grps_);
733     _gmx_sel_set_lex_input_str(scanner, str.c_str());
734     return runParser(scanner, false, -1, std::string());
735 }
736
737
738 void
739 SelectionCollection::compile()
740 {
741     if (impl_->sc_.top == NULL && requiresTopology())
742     {
743         GMX_THROW(InconsistentInputError("Selection requires topology information, but none provided"));
744     }
745     if (!impl_->bExternalGroupsSet_)
746     {
747         setIndexGroups(NULL);
748     }
749     if (impl_->debugLevel_ >= 1)
750     {
751         printTree(stderr, false);
752     }
753
754     SelectionCompiler compiler;
755     compiler.compile(this);
756
757     if (impl_->debugLevel_ >= 1)
758     {
759         std::fprintf(stderr, "\n");
760         printTree(stderr, false);
761         std::fprintf(stderr, "\n");
762         impl_->sc_.pcc.printTree(stderr);
763         std::fprintf(stderr, "\n");
764     }
765     impl_->sc_.pcc.initEvaluation();
766     if (impl_->debugLevel_ >= 1)
767     {
768         impl_->sc_.pcc.printTree(stderr);
769         std::fprintf(stderr, "\n");
770     }
771
772     // TODO: It would be nicer to associate the name of the selection option
773     // (if available) to the error message.
774     SelectionDataList::const_iterator iter;
775     for (iter = impl_->sc_.sel.begin(); iter != impl_->sc_.sel.end(); ++iter)
776     {
777         const internal::SelectionData &sel = **iter;
778         if (sel.hasFlag(efSelection_OnlyAtoms))
779         {
780             if (!sel.hasOnlyAtoms())
781             {
782                 std::string message = formatString(
783                             "Selection '%s' does not evaluate to individual atoms. "
784                             "This is not allowed in this context.",
785                             sel.selectionText());
786                 GMX_THROW(InvalidInputError(message));
787             }
788         }
789         if (sel.hasFlag(efSelection_DisallowEmpty))
790         {
791             if (sel.posCount() == 0)
792             {
793                 std::string message = formatString(
794                             "Selection '%s' never matches any atoms.",
795                             sel.selectionText());
796                 GMX_THROW(InvalidInputError(message));
797             }
798         }
799     }
800 }
801
802
803 void
804 SelectionCollection::evaluate(t_trxframe *fr, t_pbc *pbc)
805 {
806     if (fr->natoms <= impl_->maxAtomIndex_)
807     {
808         std::string message = formatString(
809                     "Trajectory has less atoms (%d) than what is required for "
810                     "evaluating the provided selections (atoms up to index %d "
811                     "are required).", fr->natoms, impl_->maxAtomIndex_ + 1);
812         GMX_THROW(InconsistentInputError(message));
813     }
814     impl_->sc_.pcc.initFrame();
815
816     SelectionEvaluator evaluator;
817     evaluator.evaluate(this, fr, pbc);
818
819     if (impl_->debugLevel_ >= 3)
820     {
821         std::fprintf(stderr, "\n");
822         printTree(stderr, true);
823     }
824 }
825
826
827 void
828 SelectionCollection::evaluateFinal(int nframes)
829 {
830     SelectionEvaluator evaluator;
831     evaluator.evaluateFinal(this, nframes);
832 }
833
834
835 void
836 SelectionCollection::printTree(FILE *fp, bool bValues) const
837 {
838     SelectionTreeElementPointer sel = impl_->sc_.root;
839     while (sel)
840     {
841         _gmx_selelem_print_tree(fp, *sel, bValues, 0);
842         sel = sel->next;
843     }
844 }
845
846
847 void
848 SelectionCollection::printXvgrInfo(FILE *out, output_env_t oenv) const
849 {
850     if (output_env_get_xvg_format(oenv) != exvgNONE)
851     {
852         const gmx_ana_selcollection_t &sc = impl_->sc_;
853         std::fprintf(out, "# Selections:\n");
854         for (int i = 0; i < sc.nvars; ++i)
855         {
856             std::fprintf(out, "#   %s\n", sc.varstrs[i]);
857         }
858         for (size_t i = 0; i < sc.sel.size(); ++i)
859         {
860             std::fprintf(out, "#   %s\n", sc.sel[i]->selectionText());
861         }
862         std::fprintf(out, "#\n");
863     }
864 }
865
866 } // namespace gmx