6b96478b87f895327060d9279e62857bd33d914f
[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,2016,2017,2018, 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 <memory>
50 #include <string>
51 #include <vector>
52
53 #include "gromacs/compat/make_unique.h"
54 #include "gromacs/onlinehelp/helpmanager.h"
55 #include "gromacs/onlinehelp/helpwritercontext.h"
56 #include "gromacs/options/basicoptions.h"
57 #include "gromacs/options/ioptionscontainer.h"
58 #include "gromacs/selection/selection.h"
59 #include "gromacs/selection/selhelp.h"
60 #include "gromacs/topology/mtop_util.h"
61 #include "gromacs/topology/topology.h"
62 #include "gromacs/trajectory/trajectoryframe.h"
63 #include "gromacs/utility/exceptions.h"
64 #include "gromacs/utility/filestream.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/smalloc.h"
67 #include "gromacs/utility/stringutil.h"
68 #include "gromacs/utility/textwriter.h"
69
70 #include "compiler.h"
71 #include "mempool.h"
72 #include "parser.h"
73 #include "poscalc.h"
74 #include "scanner.h"
75 #include "selectioncollection-impl.h"
76 #include "selelem.h"
77 #include "selmethod.h"
78 #include "symrec.h"
79
80 namespace gmx
81 {
82
83 /********************************************************************
84  * SelectionCollection::Impl
85  */
86
87 SelectionCollection::Impl::Impl()
88     : debugLevel_(0), bExternalGroupsSet_(false), grps_(nullptr)
89 {
90     sc_.nvars     = 0;
91     sc_.varstrs   = nullptr;
92     sc_.top       = nullptr;
93     gmx_ana_index_clear(&sc_.gall);
94     sc_.mempool   = nullptr;
95     sc_.symtab    = compat::make_unique<SelectionParserSymbolTable>();
96     gmx_ana_index_clear(&requiredAtoms_);
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     gmx_ana_index_deinit(&requiredAtoms_);
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]  inputStream   Stream to read from (typically the StandardInputStream).
136  * \param[in]  statusWriter  Stream to print prompts to (if NULL, no output is done).
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(TextInputStream *inputStream, TextWriter *statusWriter,
143                 std::string *line)
144 {
145     if (statusWriter != nullptr)
146     {
147         statusWriter->writeString("> ");
148     }
149     if (!inputStream->readLine(line))
150     {
151         return false;
152     }
153     while (endsWith(*line, "\\\n"))
154     {
155         line->resize(line->length() - 2);
156         if (statusWriter != nullptr)
157         {
158             statusWriter->writeString("... ");
159         }
160         std::string buffer;
161         // Return value ignored, buffer remains empty and works correctly
162         // if there is nothing to read.
163         inputStream->readLine(&buffer);
164         line->append(buffer);
165     }
166     if (endsWith(*line, "\n"))
167     {
168         line->resize(line->length() - 1);
169     }
170     else if (statusWriter != nullptr)
171     {
172         statusWriter->writeLine();
173     }
174     return true;
175 }
176
177 /*! \brief
178  * Helper function for tokenizing the input and pushing them to the parser.
179  *
180  * \param     scanner       Tokenizer data structure.
181  * \param     parserState   Parser data structure.
182  * \param[in] bInteractive  Whether to operate in interactive mode.
183  *
184  * Repeatedly reads tokens using \p scanner and pushes them to the parser with
185  * \p parserState until there is no more input, or until enough input is given
186  * (only in interactive mode).
187  */
188 int runParserLoop(yyscan_t scanner, _gmx_sel_yypstate *parserState,
189                   bool bInteractive)
190 {
191     int status    = YYPUSH_MORE;
192     do
193     {
194         YYSTYPE value;
195         YYLTYPE location;
196         int     token = _gmx_sel_yylex(&value, &location, scanner);
197         if (bInteractive && token == 0)
198         {
199             break;
200         }
201         status = _gmx_sel_yypush_parse(parserState, token, &value, &location, scanner);
202     }
203     while (status == YYPUSH_MORE);
204     _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
205     return status;
206 }
207
208 /*! \brief
209  * Print current status in response to empty line in interactive input.
210  *
211  * \param[in] writer         Writer to use for the output.
212  * \param[in] sc             Selection collection data structure.
213  * \param[in] grps           Available index groups.
214  * \param[in] firstSelection Index of first selection from this interactive
215  *     session.
216  * \param[in] maxCount       Maximum number of selections.
217  * \param[in] context        Context to print for what the selections are for.
218  * \param[in] bFirst         Whether this is the header that is printed before
219  *     any user input.
220  *
221  * Prints the available index groups and currently provided selections.
222  */
223 void printCurrentStatus(TextWriter *writer, gmx_ana_selcollection_t *sc,
224                         gmx_ana_indexgrps_t *grps, size_t firstSelection,
225                         int maxCount, const std::string &context, bool bFirst)
226 {
227     if (grps != nullptr)
228     {
229         writer->writeLine("Available static index groups:");
230         gmx_ana_indexgrps_print(writer, grps, 0);
231     }
232     writer->writeString("Specify ");
233     if (maxCount < 0)
234     {
235         writer->writeString("any number of selections");
236     }
237     else if (maxCount == 1)
238     {
239         writer->writeString("a selection");
240     }
241     else
242     {
243         writer->writeString(formatString("%d selections", maxCount));
244     }
245     writer->writeString(formatString("%s%s:\n",
246                                      context.empty() ? "" : " ", context.c_str()));
247     writer->writeString(formatString(
248                                 "(one per line, <enter> for status/groups, 'help' for help%s)\n",
249                                 maxCount < 0 ? ", Ctrl-D to end" : ""));
250     if (!bFirst && (sc->nvars > 0 || sc->sel.size() > firstSelection))
251     {
252         writer->writeLine("Currently provided selections:");
253         for (int i = 0; i < sc->nvars; ++i)
254         {
255             writer->writeString(formatString("     %s\n", sc->varstrs[i]));
256         }
257         for (size_t i = firstSelection; i < sc->sel.size(); ++i)
258         {
259             writer->writeString(formatString(
260                                         " %2d. %s\n",
261                                         static_cast<int>(i - firstSelection + 1),
262                                         sc->sel[i]->selectionText()));
263         }
264         if (maxCount > 0)
265         {
266             const int remaining
267                 = maxCount - static_cast<int>(sc->sel.size() - firstSelection);
268             writer->writeString(formatString(
269                                         "(%d more selection%s required)\n",
270                                         remaining, remaining > 1 ? "s" : ""));
271         }
272     }
273 }
274
275 /*! \brief
276  * Prints selection help in interactive selection input.
277  *
278  * \param[in] writer Writer to use for the output.
279  * \param[in] sc    Selection collection data structure.
280  * \param[in] line  Line of user input requesting help (starting with `help`).
281  *
282  * Initializes the selection help if not yet initialized, and finds the help
283  * topic based on words on the input line.
284  */
285 void printHelp(TextWriter *writer, gmx_ana_selcollection_t *sc,
286                const std::string &line)
287 {
288     if (sc->rootHelp.get() == nullptr)
289     {
290         sc->rootHelp = createSelectionHelpTopic();
291     }
292     HelpWriterContext context(writer, eHelpOutputFormat_Console);
293     HelpManager       manager(*sc->rootHelp, context);
294     try
295     {
296         std::vector<std::string>                 topic = splitString(line);
297         std::vector<std::string>::const_iterator value;
298         // First item in the list is the 'help' token.
299         for (value = topic.begin() + 1; value != topic.end(); ++value)
300         {
301             manager.enterTopic(*value);
302         }
303     }
304     catch (const InvalidInputError &ex)
305     {
306         writer->writeLine(ex.what());
307         return;
308     }
309     manager.writeCurrentTopic();
310 }
311
312 /*! \brief
313  * Helper function that runs the parser once the tokenizer has been
314  * initialized.
315  *
316  * \param[in,out] scanner       Scanner data structure.
317  * \param[in]     inputStream   Stream to use for input (currently only with
318  *      `bInteractive==true`).
319  * \param[in]     bInteractive  Whether to use a line-based reading
320  *      algorithm designed for interactive input.
321  * \param[in]     maxnr   Maximum number of selections to parse
322  *      (if -1, parse as many as provided by the user).
323  * \param[in]     context Context to print for what the selections are for.
324  * \returns       Vector of parsed selections.
325  * \throws        std::bad_alloc if out of memory.
326  * \throws        InvalidInputError if there is a parsing error.
327  *
328  * Used internally to implement parseInteractive(), parseFromFile() and
329  * parseFromString().
330  */
331 SelectionList runParser(yyscan_t scanner, TextInputStream *inputStream,
332                         bool bInteractive, int maxnr, const std::string &context)
333 {
334     std::shared_ptr<void>    scannerGuard(scanner, &_gmx_sel_free_lexer);
335     gmx_ana_selcollection_t *sc   = _gmx_sel_lexer_selcollection(scanner);
336     gmx_ana_indexgrps_t     *grps = _gmx_sel_lexer_indexgrps(scanner);
337
338     size_t                   oldCount = sc->sel.size();
339     {
340         std::shared_ptr<_gmx_sel_yypstate> parserState(
341                 _gmx_sel_yypstate_new(), &_gmx_sel_yypstate_delete);
342         if (bInteractive)
343         {
344             TextWriter *statusWriter = _gmx_sel_lexer_get_status_writer(scanner);
345             if (statusWriter != nullptr)
346             {
347                 printCurrentStatus(statusWriter, sc, grps, oldCount, maxnr, context, true);
348             }
349             std::string line;
350             int         status;
351             while (promptLine(inputStream, statusWriter, &line))
352             {
353                 if (statusWriter != nullptr)
354                 {
355                     line = stripString(line);
356                     if (line.empty())
357                     {
358                         printCurrentStatus(statusWriter, sc, grps, oldCount, maxnr, context, false);
359                         continue;
360                     }
361                     if (startsWith(line, "help")
362                         && (line[4] == 0 || (std::isspace(line[4]) != 0)))
363                     {
364                         printHelp(statusWriter, sc, line);
365                         continue;
366                     }
367                 }
368                 line.append("\n");
369                 _gmx_sel_set_lex_input_str(scanner, line.c_str());
370                 status = runParserLoop(scanner, parserState.get(), true);
371                 if (status != YYPUSH_MORE)
372                 {
373                     // TODO: Check if there is more input, and issue an
374                     // error/warning if some input was ignored.
375                     goto early_termination;
376                 }
377             }
378             {
379                 YYLTYPE location;
380                 status = _gmx_sel_yypush_parse(parserState.get(), 0, nullptr,
381                                                &location, scanner);
382             }
383             // TODO: Remove added selections from the collection if parsing failed?
384             _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
385 early_termination:
386             GMX_RELEASE_ASSERT(status == 0,
387                                "Parser errors should have resulted in an exception");
388         }
389         else
390         {
391             int status = runParserLoop(scanner, parserState.get(), false);
392             GMX_RELEASE_ASSERT(status == 0,
393                                "Parser errors should have resulted in an exception");
394         }
395     }
396     scannerGuard.reset();
397     int nr = sc->sel.size() - oldCount;
398     if (maxnr > 0 && nr != maxnr)
399     {
400         std::string message
401             = formatString("Too few selections provided; got %d, expected %d",
402                            nr, maxnr);
403         GMX_THROW(InvalidInputError(message));
404     }
405
406     SelectionList                     result;
407     SelectionDataList::const_iterator i;
408     result.reserve(nr);
409     for (i = sc->sel.begin() + oldCount; i != sc->sel.end(); ++i)
410     {
411         result.emplace_back(i->get());
412     }
413     return result;
414 }
415
416 /*! \brief
417  * Checks that index groups have valid atom indices.
418  *
419  * \param[in]    root    Root of selection tree to process.
420  * \param[in]    natoms  Maximum number of atoms that the selections are set
421  *     to evaluate.
422  * \param        errors  Object for reporting any error messages.
423  * \throws std::bad_alloc if out of memory.
424  *
425  * Recursively checks the selection tree for index groups.
426  * Each found group is checked that it only contains atom indices that match
427  * the topology/maximum number of atoms set for the selection collection.
428  * Any issues are reported to \p errors.
429  */
430 void checkExternalGroups(const SelectionTreeElementPointer &root,
431                          int                                natoms,
432                          ExceptionInitializer              *errors)
433 {
434     if (root->type == SEL_CONST && root->v.type == GROUP_VALUE)
435     {
436         try
437         {
438             root->checkIndexGroup(natoms);
439         }
440         catch (const UserInputError &)
441         {
442             errors->addCurrentExceptionAsNested();
443         }
444     }
445
446     SelectionTreeElementPointer child = root->child;
447     while (child)
448     {
449         checkExternalGroups(child, natoms, errors);
450         child = child->next;
451     }
452 }
453
454 //! Checks whether the given topology properties are available.
455 void checkTopologyProperties(const gmx_mtop_t                  *top,
456                              const SelectionTopologyProperties &props)
457 {
458     if (top == nullptr)
459     {
460         if (props.hasAny())
461         {
462             GMX_THROW(InconsistentInputError("Selection requires topology information, but none provided"));
463         }
464         return;
465     }
466     if (props.needsMasses && !gmx_mtop_has_masses(top))
467     {
468         GMX_THROW(InconsistentInputError("Selection requires mass information, but it is not available in the topology"));
469     }
470 }
471
472 }   // namespace
473
474
475 void SelectionCollection::Impl::resolveExternalGroups(
476         const SelectionTreeElementPointer &root,
477         ExceptionInitializer              *errors)
478 {
479
480     if (root->type == SEL_GROUPREF)
481     {
482         try
483         {
484             root->resolveIndexGroupReference(grps_, sc_.gall.isize);
485         }
486         catch (const UserInputError &)
487         {
488             errors->addCurrentExceptionAsNested();
489         }
490     }
491
492     SelectionTreeElementPointer child = root->child;
493     while (child)
494     {
495         resolveExternalGroups(child, errors);
496         root->flags |= (child->flags & SEL_UNSORTED);
497         child        = child->next;
498     }
499 }
500
501
502 bool SelectionCollection::Impl::areForcesRequested() const
503 {
504     for (const auto &sel : sc_.sel)
505     {
506         if (sel->hasFlag(gmx::efSelection_EvaluateForces))
507         {
508             return true;
509         }
510     }
511     return false;
512 }
513
514
515 SelectionTopologyProperties
516 SelectionCollection::Impl::requiredTopologyPropertiesForPositionType(
517         const std::string &post, bool forces) const
518 {
519     SelectionTopologyProperties props;
520     if (!post.empty())
521     {
522         switch (PositionCalculationCollection::requiredTopologyInfoForType(post.c_str(), forces))
523         {
524             case PositionCalculationCollection::RequiredTopologyInfo::None:
525                 break;
526             case PositionCalculationCollection::RequiredTopologyInfo::Topology:
527                 props.merge(SelectionTopologyProperties::topology());
528                 break;
529             case PositionCalculationCollection::RequiredTopologyInfo::TopologyAndMasses:
530                 props.merge(SelectionTopologyProperties::masses());
531                 break;
532         }
533     }
534     return props;
535 }
536
537
538 /********************************************************************
539  * SelectionCollection
540  */
541
542 SelectionCollection::SelectionCollection()
543     : impl_(new Impl)
544 {
545 }
546
547
548 SelectionCollection::~SelectionCollection()
549 {
550 }
551
552
553 void
554 SelectionCollection::initOptions(IOptionsContainer   *options,
555                                  SelectionTypeOption  selectionTypeOption)
556 {
557     const char * const debug_levels[]
558         = { "no", "basic", "compile", "eval", "full" };
559
560     const char *const *postypes = PositionCalculationCollection::typeEnumValues;
561     options->addOption(StringOption("selrpos")
562                            .enumValueFromNullTerminatedArray(postypes)
563                            .store(&impl_->rpost_).defaultValue(postypes[0])
564                            .description("Selection reference positions"));
565     if (selectionTypeOption == IncludeSelectionTypeOption)
566     {
567         options->addOption(StringOption("seltype")
568                                .enumValueFromNullTerminatedArray(postypes)
569                                .store(&impl_->spost_).defaultValue(postypes[0])
570                                .description("Default selection output positions"));
571     }
572     else
573     {
574         impl_->spost_ = postypes[0];
575     }
576     GMX_RELEASE_ASSERT(impl_->debugLevel_ >= 0 && impl_->debugLevel_ <= 4,
577                        "Debug level out of range");
578     options->addOption(EnumIntOption("seldebug").hidden(impl_->debugLevel_ == 0)
579                            .enumValue(debug_levels).store(&impl_->debugLevel_)
580                            .description("Print out selection trees for debugging"));
581 }
582
583
584 void
585 SelectionCollection::setReferencePosType(const char *type)
586 {
587     GMX_RELEASE_ASSERT(type != nullptr, "Cannot assign NULL position type");
588     // Check that the type is valid, throw if it is not.
589     e_poscalc_t  dummytype;
590     int          dummyflags;
591     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
592     impl_->rpost_ = type;
593 }
594
595
596 void
597 SelectionCollection::setOutputPosType(const char *type)
598 {
599     GMX_RELEASE_ASSERT(type != nullptr, "Cannot assign NULL position type");
600     // Check that the type is valid, throw if it is not.
601     e_poscalc_t  dummytype;
602     int          dummyflags;
603     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
604     impl_->spost_ = type;
605 }
606
607
608 void
609 SelectionCollection::setDebugLevel(int debugLevel)
610 {
611     impl_->debugLevel_ = debugLevel;
612 }
613
614
615 void
616 SelectionCollection::setTopology(gmx_mtop_t *top, int natoms)
617 {
618     GMX_RELEASE_ASSERT(natoms > 0 || top != nullptr,
619                        "The number of atoms must be given if there is no topology");
620     checkTopologyProperties(top, requiredTopologyProperties());
621     // Get the number of atoms from the topology if it is not given.
622     if (natoms <= 0)
623     {
624         natoms = top->natoms;
625     }
626     if (impl_->bExternalGroupsSet_)
627     {
628         ExceptionInitializer        errors("Invalid index group references encountered");
629         SelectionTreeElementPointer root = impl_->sc_.root;
630         while (root)
631         {
632             checkExternalGroups(root, natoms, &errors);
633             root = root->next;
634         }
635         if (errors.hasNestedExceptions())
636         {
637             GMX_THROW(InconsistentInputError(errors));
638         }
639     }
640     gmx_ana_selcollection_t *sc = &impl_->sc_;
641     // Do this first, as it allocates memory, while the others don't throw.
642     gmx_ana_index_init_simple(&sc->gall, natoms);
643     sc->top = top;
644     sc->pcc.setTopology(top);
645 }
646
647
648 void
649 SelectionCollection::setIndexGroups(gmx_ana_indexgrps_t *grps)
650 {
651     GMX_RELEASE_ASSERT(grps == nullptr || !impl_->bExternalGroupsSet_,
652                        "Can only set external groups once or clear them afterwards");
653     impl_->grps_               = grps;
654     impl_->bExternalGroupsSet_ = true;
655
656     ExceptionInitializer        errors("Invalid index group reference(s)");
657     SelectionTreeElementPointer root = impl_->sc_.root;
658     while (root)
659     {
660         impl_->resolveExternalGroups(root, &errors);
661         root->checkUnsortedAtoms(true, &errors);
662         root = root->next;
663     }
664     if (errors.hasNestedExceptions())
665     {
666         GMX_THROW(InconsistentInputError(errors));
667     }
668     for (size_t i = 0; i < impl_->sc_.sel.size(); ++i)
669     {
670         impl_->sc_.sel[i]->refreshName();
671     }
672 }
673
674 SelectionTopologyProperties
675 SelectionCollection::requiredTopologyProperties() const
676 {
677     SelectionTopologyProperties props;
678
679     // These should not throw, because has been checked earlier.
680     props.merge(impl_->requiredTopologyPropertiesForPositionType(impl_->rpost_, false));
681     const bool forcesRequested = impl_->areForcesRequested();
682     props.merge(impl_->requiredTopologyPropertiesForPositionType(impl_->spost_,
683                                                                  forcesRequested));
684
685     SelectionTreeElementPointer sel = impl_->sc_.root;
686     while (sel && !props.hasAll())
687     {
688         props.merge(sel->requiredTopologyProperties());
689         sel = sel->next;
690     }
691     return props;
692 }
693
694
695 bool
696 SelectionCollection::requiresIndexGroups() const
697 {
698     SelectionTreeElementPointer sel = impl_->sc_.root;
699     while (sel)
700     {
701         if (sel->requiresIndexGroups())
702         {
703             return true;
704         }
705         sel = sel->next;
706     }
707     return false;
708 }
709
710
711 SelectionList
712 SelectionCollection::parseFromStdin(int count, bool bInteractive,
713                                     const std::string &context)
714 {
715     return parseInteractive(count, &StandardInputStream::instance(),
716                             bInteractive ? &TextOutputFile::standardError() : nullptr,
717                             context);
718 }
719
720 namespace
721 {
722
723 //! Helper function to initialize status writer for interactive selection parsing.
724 std::unique_ptr<TextWriter> initStatusWriter(TextOutputStream *statusStream)
725 {
726     std::unique_ptr<TextWriter> statusWriter;
727     if (statusStream != nullptr)
728     {
729         statusWriter = compat::make_unique<TextWriter>(statusStream);
730         statusWriter->wrapperSettings().setLineLength(78);
731     }
732     return statusWriter;
733 }
734
735 }   // namespace
736
737 SelectionList
738 SelectionCollection::parseInteractive(int                count,
739                                       TextInputStream   *inputStream,
740                                       TextOutputStream  *statusStream,
741                                       const std::string &context)
742 {
743     yyscan_t scanner;
744
745     const std::unique_ptr<TextWriter> statusWriter(
746             initStatusWriter(statusStream));
747     _gmx_sel_init_lexer(&scanner, &impl_->sc_, statusWriter.get(),
748                         count, impl_->bExternalGroupsSet_, impl_->grps_);
749     return runParser(scanner, inputStream, true, count, context);
750 }
751
752
753 SelectionList
754 SelectionCollection::parseFromFile(const std::string &filename)
755 {
756
757     try
758     {
759         yyscan_t      scanner;
760         TextInputFile file(filename);
761         // TODO: Exception-safe way of using the lexer.
762         _gmx_sel_init_lexer(&scanner, &impl_->sc_, nullptr, -1,
763                             impl_->bExternalGroupsSet_,
764                             impl_->grps_);
765         _gmx_sel_set_lex_input_file(scanner, file.handle());
766         return runParser(scanner, nullptr, false, -1, std::string());
767     }
768     catch (GromacsException &ex)
769     {
770         ex.prependContext(formatString(
771                                   "Error in parsing selections from file '%s'",
772                                   filename.c_str()));
773         throw;
774     }
775 }
776
777
778 SelectionList
779 SelectionCollection::parseFromString(const std::string &str)
780 {
781     yyscan_t scanner;
782
783     _gmx_sel_init_lexer(&scanner, &impl_->sc_, nullptr, -1,
784                         impl_->bExternalGroupsSet_,
785                         impl_->grps_);
786     _gmx_sel_set_lex_input_str(scanner, str.c_str());
787     return runParser(scanner, nullptr, false, -1, std::string());
788 }
789
790
791 void
792 SelectionCollection::compile()
793 {
794     checkTopologyProperties(impl_->sc_.top, requiredTopologyProperties());
795     if (!impl_->bExternalGroupsSet_)
796     {
797         setIndexGroups(nullptr);
798     }
799     if (impl_->debugLevel_ >= 1)
800     {
801         printTree(stderr, false);
802     }
803
804     SelectionCompiler compiler;
805     compiler.compile(this);
806
807     if (impl_->debugLevel_ >= 1)
808     {
809         std::fprintf(stderr, "\n");
810         printTree(stderr, false);
811         std::fprintf(stderr, "\n");
812         impl_->sc_.pcc.printTree(stderr);
813         std::fprintf(stderr, "\n");
814     }
815     impl_->sc_.pcc.initEvaluation();
816     if (impl_->debugLevel_ >= 1)
817     {
818         impl_->sc_.pcc.printTree(stderr);
819         std::fprintf(stderr, "\n");
820     }
821
822     // TODO: It would be nicer to associate the name of the selection option
823     // (if available) to the error message.
824     SelectionDataList::const_iterator iter;
825     for (iter = impl_->sc_.sel.begin(); iter != impl_->sc_.sel.end(); ++iter)
826     {
827         const internal::SelectionData &sel = **iter;
828         if (sel.hasFlag(efSelection_OnlyAtoms))
829         {
830             if (!sel.hasOnlyAtoms())
831             {
832                 std::string message = formatString(
833                             "Selection '%s' does not evaluate to individual atoms. "
834                             "This is not allowed in this context.",
835                             sel.selectionText());
836                 GMX_THROW(InvalidInputError(message));
837             }
838             if (sel.hasFlag(efSelection_OnlySorted))
839             {
840                 if (!sel.hasSortedAtomIndices())
841                 {
842                     const std::string message = formatString(
843                                 "Selection '%s' does not evaluate to atoms in an "
844                                 "ascending (sorted) order. "
845                                 "This is not allowed in this context.",
846                                 sel.selectionText());
847                     GMX_THROW(InvalidInputError(message));
848                 }
849             }
850         }
851         if (sel.hasFlag(efSelection_DisallowEmpty))
852         {
853             if (sel.posCount() == 0)
854             {
855                 std::string message = formatString(
856                             "Selection '%s' never matches any atoms.",
857                             sel.selectionText());
858                 GMX_THROW(InvalidInputError(message));
859             }
860         }
861     }
862     impl_->rpost_.clear();
863     impl_->spost_.clear();
864 }
865
866
867 void
868 SelectionCollection::evaluate(t_trxframe *fr, t_pbc *pbc)
869 {
870     checkTopologyProperties(impl_->sc_.top, requiredTopologyProperties());
871     if (fr->bIndex)
872     {
873         gmx_ana_index_t g;
874         gmx_ana_index_set(&g, fr->natoms, fr->index, 0);
875         GMX_RELEASE_ASSERT(gmx_ana_index_check_sorted(&g),
876                            "Only trajectories with atoms in ascending order "
877                            "are currently supported");
878         if (!gmx_ana_index_contains(&g, &impl_->requiredAtoms_))
879         {
880             const std::string message = formatString(
881                         "Trajectory does not contain all atoms required for "
882                         "evaluating the provided selections.");
883             GMX_THROW(InconsistentInputError(message));
884         }
885     }
886     else
887     {
888         const int maxAtomIndex = gmx_ana_index_get_max_index(&impl_->requiredAtoms_);
889         if (fr->natoms <= maxAtomIndex)
890         {
891             const std::string message = formatString(
892                         "Trajectory has less atoms (%d) than what is required for "
893                         "evaluating the provided selections (atoms up to index %d "
894                         "are required).", fr->natoms, maxAtomIndex + 1);
895             GMX_THROW(InconsistentInputError(message));
896         }
897     }
898     impl_->sc_.pcc.initFrame(fr);
899
900     SelectionEvaluator evaluator;
901     evaluator.evaluate(this, fr, pbc);
902
903     if (impl_->debugLevel_ >= 3)
904     {
905         std::fprintf(stderr, "\n");
906         printTree(stderr, true);
907     }
908 }
909
910
911 void
912 SelectionCollection::evaluateFinal(int nframes)
913 {
914     SelectionEvaluator evaluator;
915     evaluator.evaluateFinal(this, nframes);
916 }
917
918
919 void
920 SelectionCollection::printTree(FILE *fp, bool bValues) const
921 {
922     SelectionTreeElementPointer sel = impl_->sc_.root;
923     while (sel)
924     {
925         _gmx_selelem_print_tree(fp, *sel, bValues, 0);
926         sel = sel->next;
927     }
928 }
929
930
931 void
932 SelectionCollection::printXvgrInfo(FILE *out) const
933 {
934     const gmx_ana_selcollection_t &sc = impl_->sc_;
935     std::fprintf(out, "# Selections:\n");
936     for (int i = 0; i < sc.nvars; ++i)
937     {
938         std::fprintf(out, "#   %s\n", sc.varstrs[i]);
939     }
940     for (size_t i = 0; i < sc.sel.size(); ++i)
941     {
942         std::fprintf(out, "#   %s\n", sc.sel[i]->selectionText());
943     }
944     std::fprintf(out, "#\n");
945 }
946
947 } // namespace gmx