d124cf2800bcc15109af43e2dcc3f4caa375b1bf
[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,2015, 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/scoped_ptr.hpp>
53 #include <boost/shared_ptr.hpp>
54
55 #include "gromacs/fileio/trx.h"
56 #include "gromacs/legacyheaders/oenv.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/selection/selhelp.h"
63 #include "gromacs/topology/topology.h"
64 #include "gromacs/utility/exceptions.h"
65 #include "gromacs/utility/filestream.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/smalloc.h"
68 #include "gromacs/utility/stringutil.h"
69 #include "gromacs/utility/textwriter.h"
70
71 #include "compiler.h"
72 #include "mempool.h"
73 #include "parser.h"
74 #include "poscalc.h"
75 #include "scanner.h"
76 #include "selectioncollection-impl.h"
77 #include "selelem.h"
78 #include "selmethod.h"
79 #include "symrec.h"
80
81 namespace gmx
82 {
83
84 /********************************************************************
85  * SelectionCollection::Impl
86  */
87
88 SelectionCollection::Impl::Impl()
89     : maxAtomIndex_(0), debugLevel_(0), bExternalGroupsSet_(false), grps_(NULL)
90 {
91     sc_.nvars     = 0;
92     sc_.varstrs   = NULL;
93     sc_.top       = NULL;
94     gmx_ana_index_clear(&sc_.gall);
95     sc_.mempool   = NULL;
96     sc_.symtab.reset(new SelectionParserSymbolTable);
97     gmx_ana_selmethod_register_defaults(sc_.symtab.get());
98 }
99
100
101 SelectionCollection::Impl::~Impl()
102 {
103     clearSymbolTable();
104     // The tree must be freed before the SelectionData objects, since the
105     // tree may hold references to the position data in SelectionData.
106     sc_.root.reset();
107     sc_.sel.clear();
108     for (int i = 0; i < sc_.nvars; ++i)
109     {
110         sfree(sc_.varstrs[i]);
111     }
112     sfree(sc_.varstrs);
113     gmx_ana_index_deinit(&sc_.gall);
114     if (sc_.mempool)
115     {
116         _gmx_sel_mempool_destroy(sc_.mempool);
117     }
118 }
119
120
121 void
122 SelectionCollection::Impl::clearSymbolTable()
123 {
124     sc_.symtab.reset();
125 }
126
127
128 namespace
129 {
130
131 /*! \brief
132  * Reads a single selection line from stdin.
133  *
134  * \param[in]  inputStream   Stream to read from (typically the StandardInputStream).
135  * \param[in]  statusWriter  Stream to print prompts to (if NULL, no output is done).
136  * \param[out] line          The read line in stored here.
137  * \returns true if something was read, false if at end of input.
138  *
139  * Handles line continuation, reading also the continuing line(s) in one call.
140  */
141 bool promptLine(TextInputStream *inputStream, TextWriter *statusWriter,
142                 std::string *line)
143 {
144     if (statusWriter != NULL)
145     {
146         statusWriter->writeString("> ");
147     }
148     if (!inputStream->readLine(line))
149     {
150         return false;
151     }
152     while (endsWith(*line, "\\\n"))
153     {
154         line->resize(line->length() - 2);
155         if (statusWriter != NULL)
156         {
157             statusWriter->writeString("... ");
158         }
159         std::string buffer;
160         // Return value ignored, buffer remains empty and works correctly
161         // if there is nothing to read.
162         inputStream->readLine(&buffer);
163         line->append(buffer);
164     }
165     if (endsWith(*line, "\n"))
166     {
167         line->resize(line->length() - 1);
168     }
169     else if (statusWriter != NULL)
170     {
171         statusWriter->writeLine();
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     do
192     {
193         YYSTYPE value;
194         YYLTYPE location;
195         int     token = _gmx_sel_yylex(&value, &location, scanner);
196         if (bInteractive && token == 0)
197         {
198             break;
199         }
200         status = _gmx_sel_yypush_parse(parserState, token, &value, &location, scanner);
201     }
202     while (status == YYPUSH_MORE);
203     _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
204     return status;
205 }
206
207 /*! \brief
208  * Print current status in response to empty line in interactive input.
209  *
210  * \param[in] writer         Writer to use for the output.
211  * \param[in] sc             Selection collection data structure.
212  * \param[in] grps           Available index groups.
213  * \param[in] firstSelection Index of first selection from this interactive
214  *     session.
215  * \param[in] maxCount       Maximum number of selections.
216  * \param[in] context        Context to print for what the selections are for.
217  * \param[in] bFirst         Whether this is the header that is printed before
218  *     any user input.
219  *
220  * Prints the available index groups and currently provided selections.
221  */
222 void printCurrentStatus(TextWriter *writer, gmx_ana_selcollection_t *sc,
223                         gmx_ana_indexgrps_t *grps, size_t firstSelection,
224                         int maxCount, const std::string &context, bool bFirst)
225 {
226     if (grps != NULL)
227     {
228         writer->writeLine("Available static index groups:");
229         gmx_ana_indexgrps_print(writer, grps, 0);
230     }
231     writer->writeString("Specify ");
232     if (maxCount < 0)
233     {
234         writer->writeString("any number of selections");
235     }
236     else if (maxCount == 1)
237     {
238         writer->writeString("a selection");
239     }
240     else
241     {
242         writer->writeString(formatString("%d selections", maxCount));
243     }
244     writer->writeString(formatString("%s%s:\n",
245                                      context.empty() ? "" : " ", context.c_str()));
246     writer->writeString(formatString(
247                                 "(one per line, <enter> for status/groups, 'help' for help%s)\n",
248                                 maxCount < 0 ? ", Ctrl-D to end" : ""));
249     if (!bFirst && (sc->nvars > 0 || sc->sel.size() > firstSelection))
250     {
251         writer->writeLine("Currently provided selections:");
252         for (int i = 0; i < sc->nvars; ++i)
253         {
254             writer->writeString(formatString("     %s\n", sc->varstrs[i]));
255         }
256         for (size_t i = firstSelection; i < sc->sel.size(); ++i)
257         {
258             writer->writeString(formatString(
259                                         " %2d. %s\n",
260                                         static_cast<int>(i - firstSelection + 1),
261                                         sc->sel[i]->selectionText()));
262         }
263         if (maxCount > 0)
264         {
265             const int remaining
266                 = maxCount - static_cast<int>(sc->sel.size() - firstSelection);
267             writer->writeString(formatString(
268                                         "(%d more selection%s required)\n",
269                                         remaining, remaining > 1 ? "s" : ""));
270         }
271     }
272 }
273
274 /*! \brief
275  * Prints selection help in interactive selection input.
276  *
277  * \param[in] writer Writer to use for the output.
278  * \param[in] sc    Selection collection data structure.
279  * \param[in] line  Line of user input requesting help (starting with `help`).
280  *
281  * Initializes the selection help if not yet initialized, and finds the help
282  * topic based on words on the input line.
283  */
284 void printHelp(TextWriter *writer, gmx_ana_selcollection_t *sc,
285                const std::string &line)
286 {
287     if (sc->rootHelp.get() == NULL)
288     {
289         sc->rootHelp = createSelectionHelpTopic();
290     }
291     HelpWriterContext context(&writer->stream(), eHelpOutputFormat_Console);
292     HelpManager       manager(*sc->rootHelp, context);
293     try
294     {
295         std::vector<std::string>                 topic = splitString(line);
296         std::vector<std::string>::const_iterator value;
297         // First item in the list is the 'help' token.
298         for (value = topic.begin() + 1; value != topic.end(); ++value)
299         {
300             manager.enterTopic(*value);
301         }
302     }
303     catch (const InvalidInputError &ex)
304     {
305         writer->writeLine(ex.what());
306         return;
307     }
308     manager.writeCurrentTopic();
309 }
310
311 /*! \brief
312  * Helper function that runs the parser once the tokenizer has been
313  * initialized.
314  *
315  * \param[in,out] scanner       Scanner data structure.
316  * \param[in]     inputStream   Stream to use for input (currently only with
317  *      `bInteractive==true`).
318  * \param[in]     bInteractive  Whether to use a line-based reading
319  *      algorithm designed for interactive input.
320  * \param[in]     maxnr   Maximum number of selections to parse
321  *      (if -1, parse as many as provided by the user).
322  * \param[in]     context Context to print for what the selections are for.
323  * \returns       Vector of parsed selections.
324  * \throws        std::bad_alloc if out of memory.
325  * \throws        InvalidInputError if there is a parsing error.
326  *
327  * Used internally to implement parseInteractive(), parseFromFile() and
328  * parseFromString().
329  */
330 SelectionList runParser(yyscan_t scanner, TextInputStream *inputStream,
331                         bool bInteractive, int maxnr, const std::string &context)
332 {
333     boost::shared_ptr<void>  scannerGuard(scanner, &_gmx_sel_free_lexer);
334     gmx_ana_selcollection_t *sc   = _gmx_sel_lexer_selcollection(scanner);
335     gmx_ana_indexgrps_t     *grps = _gmx_sel_lexer_indexgrps(scanner);
336
337     size_t                   oldCount = sc->sel.size();
338     {
339         boost::shared_ptr<_gmx_sel_yypstate> parserState(
340                 _gmx_sel_yypstate_new(), &_gmx_sel_yypstate_delete);
341         if (bInteractive)
342         {
343             TextWriter *statusWriter = _gmx_sel_lexer_get_status_writer(scanner);
344             if (statusWriter != NULL)
345             {
346                 printCurrentStatus(statusWriter, sc, grps, oldCount, maxnr, context, true);
347             }
348             std::string line;
349             int         status;
350             while (promptLine(inputStream, statusWriter, &line))
351             {
352                 if (statusWriter != NULL)
353                 {
354                     line = stripString(line);
355                     if (line.empty())
356                     {
357                         printCurrentStatus(statusWriter, sc, grps, oldCount, maxnr, context, false);
358                         continue;
359                     }
360                     if (startsWith(line, "help")
361                         && (line[4] == 0 || std::isspace(line[4])))
362                     {
363                         printHelp(statusWriter, sc, line);
364                         continue;
365                     }
366                 }
367                 line.append("\n");
368                 _gmx_sel_set_lex_input_str(scanner, line.c_str());
369                 status = runParserLoop(scanner, parserState.get(), true);
370                 if (status != YYPUSH_MORE)
371                 {
372                     // TODO: Check if there is more input, and issue an
373                     // error/warning if some input was ignored.
374                     goto early_termination;
375                 }
376             }
377             {
378                 YYLTYPE location;
379                 status = _gmx_sel_yypush_parse(parserState.get(), 0, NULL,
380                                                &location, scanner);
381             }
382             // TODO: Remove added selections from the collection if parsing failed?
383             _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
384 early_termination:
385             GMX_RELEASE_ASSERT(status == 0,
386                                "Parser errors should have resulted in an exception");
387         }
388         else
389         {
390             int status = runParserLoop(scanner, parserState.get(), false);
391             GMX_RELEASE_ASSERT(status == 0,
392                                "Parser errors should have resulted in an exception");
393         }
394     }
395     scannerGuard.reset();
396     int nr = sc->sel.size() - oldCount;
397     if (maxnr > 0 && nr != maxnr)
398     {
399         std::string message
400             = formatString("Too few selections provided; got %d, expected %d",
401                            nr, maxnr);
402         GMX_THROW(InvalidInputError(message));
403     }
404
405     SelectionList                     result;
406     SelectionDataList::const_iterator i;
407     result.reserve(nr);
408     for (i = sc->sel.begin() + oldCount; i != sc->sel.end(); ++i)
409     {
410         result.push_back(Selection(i->get()));
411     }
412     return result;
413 }
414
415 /*! \brief
416  * Checks that index groups have valid atom indices.
417  *
418  * \param[in]    root    Root of selection tree to process.
419  * \param[in]    natoms  Maximum number of atoms that the selections are set
420  *     to evaluate.
421  * \param        errors  Object for reporting any error messages.
422  * \throws std::bad_alloc if out of memory.
423  *
424  * Recursively checks the selection tree for index groups.
425  * Each found group is checked that it only contains atom indices that match
426  * the topology/maximum number of atoms set for the selection collection.
427  * Any issues are reported to \p errors.
428  */
429 void checkExternalGroups(const SelectionTreeElementPointer &root,
430                          int                                natoms,
431                          ExceptionInitializer              *errors)
432 {
433     if (root->type == SEL_CONST && root->v.type == GROUP_VALUE)
434     {
435         try
436         {
437             root->checkIndexGroup(natoms);
438         }
439         catch (const UserInputError &)
440         {
441             errors->addCurrentExceptionAsNested();
442         }
443     }
444
445     SelectionTreeElementPointer child = root->child;
446     while (child)
447     {
448         checkExternalGroups(child, natoms, errors);
449         child = child->next;
450     }
451 }
452
453 }   // namespace
454
455
456 void SelectionCollection::Impl::resolveExternalGroups(
457         const SelectionTreeElementPointer &root,
458         ExceptionInitializer              *errors)
459 {
460
461     if (root->type == SEL_GROUPREF)
462     {
463         try
464         {
465             root->resolveIndexGroupReference(grps_, sc_.gall.isize);
466         }
467         catch (const UserInputError &)
468         {
469             errors->addCurrentExceptionAsNested();
470         }
471     }
472
473     SelectionTreeElementPointer child = root->child;
474     while (child)
475     {
476         resolveExternalGroups(child, errors);
477         root->flags |= (child->flags & SEL_UNSORTED);
478         child        = child->next;
479     }
480 }
481
482
483 /********************************************************************
484  * SelectionCollection
485  */
486
487 SelectionCollection::SelectionCollection()
488     : impl_(new Impl)
489 {
490 }
491
492
493 SelectionCollection::~SelectionCollection()
494 {
495 }
496
497
498 void
499 SelectionCollection::initOptions(Options *options)
500 {
501     const char * const debug_levels[]
502         = { "no", "basic", "compile", "eval", "full" };
503
504     bool bAllowNonAtomOutput = false;
505     SelectionDataList::const_iterator iter;
506     for (iter = impl_->sc_.sel.begin(); iter != impl_->sc_.sel.end(); ++iter)
507     {
508         const internal::SelectionData &sel = **iter;
509         if (!sel.hasFlag(efSelection_OnlyAtoms))
510         {
511             bAllowNonAtomOutput = true;
512         }
513     }
514
515     const char *const *postypes = PositionCalculationCollection::typeEnumValues;
516     options->addOption(StringOption("selrpos")
517                            .enumValueFromNullTerminatedArray(postypes)
518                            .store(&impl_->rpost_).defaultValue(postypes[0])
519                            .description("Selection reference positions"));
520     if (bAllowNonAtomOutput)
521     {
522         options->addOption(StringOption("seltype")
523                                .enumValueFromNullTerminatedArray(postypes)
524                                .store(&impl_->spost_).defaultValue(postypes[0])
525                                .description("Default selection output positions"));
526     }
527     else
528     {
529         impl_->spost_ = postypes[0];
530     }
531     GMX_RELEASE_ASSERT(impl_->debugLevel_ >= 0 && impl_->debugLevel_ <= 4,
532                        "Debug level out of range");
533     options->addOption(StringOption("seldebug").hidden(impl_->debugLevel_ == 0)
534                            .enumValue(debug_levels)
535                            .defaultValue(debug_levels[impl_->debugLevel_])
536                            .storeEnumIndex(&impl_->debugLevel_)
537                            .description("Print out selection trees for debugging"));
538 }
539
540
541 void
542 SelectionCollection::setReferencePosType(const char *type)
543 {
544     GMX_RELEASE_ASSERT(type != NULL, "Cannot assign NULL position type");
545     // Check that the type is valid, throw if it is not.
546     e_poscalc_t  dummytype;
547     int          dummyflags;
548     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
549     impl_->rpost_ = type;
550 }
551
552
553 void
554 SelectionCollection::setOutputPosType(const char *type)
555 {
556     GMX_RELEASE_ASSERT(type != NULL, "Cannot assign NULL position type");
557     // Check that the type is valid, throw if it is not.
558     e_poscalc_t  dummytype;
559     int          dummyflags;
560     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
561     impl_->spost_ = type;
562 }
563
564
565 void
566 SelectionCollection::setDebugLevel(int debugLevel)
567 {
568     impl_->debugLevel_ = debugLevel;
569 }
570
571
572 void
573 SelectionCollection::setTopology(t_topology *top, int natoms)
574 {
575     GMX_RELEASE_ASSERT(natoms > 0 || top != NULL,
576                        "The number of atoms must be given if there is no topology");
577     // Get the number of atoms from the topology if it is not given.
578     if (natoms <= 0)
579     {
580         natoms = top->atoms.nr;
581     }
582     if (impl_->bExternalGroupsSet_)
583     {
584         ExceptionInitializer        errors("Invalid index group references encountered");
585         SelectionTreeElementPointer root = impl_->sc_.root;
586         while (root)
587         {
588             checkExternalGroups(root, natoms, &errors);
589             root = root->next;
590         }
591         if (errors.hasNestedExceptions())
592         {
593             GMX_THROW(InconsistentInputError(errors));
594         }
595     }
596     gmx_ana_selcollection_t *sc = &impl_->sc_;
597     // Do this first, as it allocates memory, while the others don't throw.
598     gmx_ana_index_init_simple(&sc->gall, natoms);
599     sc->pcc.setTopology(top);
600     sc->top = top;
601 }
602
603
604 void
605 SelectionCollection::setIndexGroups(gmx_ana_indexgrps_t *grps)
606 {
607     GMX_RELEASE_ASSERT(grps == NULL || !impl_->bExternalGroupsSet_,
608                        "Can only set external groups once or clear them afterwards");
609     impl_->grps_               = grps;
610     impl_->bExternalGroupsSet_ = true;
611
612     ExceptionInitializer        errors("Invalid index group reference(s)");
613     SelectionTreeElementPointer root = impl_->sc_.root;
614     while (root)
615     {
616         impl_->resolveExternalGroups(root, &errors);
617         root->checkUnsortedAtoms(true, &errors);
618         root = root->next;
619     }
620     if (errors.hasNestedExceptions())
621     {
622         GMX_THROW(InconsistentInputError(errors));
623     }
624     for (size_t i = 0; i < impl_->sc_.sel.size(); ++i)
625     {
626         impl_->sc_.sel[i]->refreshName();
627     }
628 }
629
630
631 bool
632 SelectionCollection::requiresTopology() const
633 {
634     e_poscalc_t  type;
635     int          flags;
636
637     if (!impl_->rpost_.empty())
638     {
639         flags = 0;
640         // Should not throw, because has been checked earlier.
641         PositionCalculationCollection::typeFromEnum(impl_->rpost_.c_str(),
642                                                     &type, &flags);
643         if (type != POS_ATOM)
644         {
645             return true;
646         }
647     }
648     if (!impl_->spost_.empty())
649     {
650         flags = 0;
651         // Should not throw, because has been checked earlier.
652         PositionCalculationCollection::typeFromEnum(impl_->spost_.c_str(),
653                                                     &type, &flags);
654         if (type != POS_ATOM)
655         {
656             return true;
657         }
658     }
659
660     SelectionTreeElementPointer sel = impl_->sc_.root;
661     while (sel)
662     {
663         if (_gmx_selelem_requires_top(*sel))
664         {
665             return true;
666         }
667         sel = sel->next;
668     }
669     return false;
670 }
671
672 SelectionList
673 SelectionCollection::parseFromStdin(int count, bool bInteractive,
674                                     const std::string &context)
675 {
676     return parseInteractive(count, &StandardInputStream::instance(),
677                             bInteractive ? &TextOutputFile::standardError() : NULL,
678                             context);
679 }
680
681 SelectionList
682 SelectionCollection::parseInteractive(int                count,
683                                       TextInputStream   *inputStream,
684                                       TextOutputStream  *statusStream,
685                                       const std::string &context)
686 {
687     yyscan_t scanner;
688
689     boost::scoped_ptr<TextWriter> statusWriter;
690     if (statusStream != NULL)
691     {
692         statusWriter.reset(new TextWriter(statusStream));
693     }
694     _gmx_sel_init_lexer(&scanner, &impl_->sc_, statusWriter.get(),
695                         count, impl_->bExternalGroupsSet_, impl_->grps_);
696     return runParser(scanner, inputStream, true, count, context);
697 }
698
699
700 SelectionList
701 SelectionCollection::parseFromFile(const std::string &filename)
702 {
703
704     try
705     {
706         yyscan_t      scanner;
707         TextInputFile file(filename);
708         // TODO: Exception-safe way of using the lexer.
709         _gmx_sel_init_lexer(&scanner, &impl_->sc_, NULL, -1,
710                             impl_->bExternalGroupsSet_,
711                             impl_->grps_);
712         _gmx_sel_set_lex_input_file(scanner, file.handle());
713         return runParser(scanner, NULL, 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_, NULL, -1,
731                         impl_->bExternalGroupsSet_,
732                         impl_->grps_);
733     _gmx_sel_set_lex_input_str(scanner, str.c_str());
734     return runParser(scanner, NULL, 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