src/gromacs/selection/selectioncollection.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010-2018, The GROMACS development team.
   5  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::SelectionCollection.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_selection
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "selectioncollection.h"
  46
  47 #include <cctype>
  48 #include <cstdio>
  49
  50 #include <memory>
  51 #include <string>
  52 #include <vector>
  53
  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_(DebugLevel::None),
  89     bExternalGroupsSet_(false),
  90     grps_(nullptr)
  91 {
  92     sc_.nvars   = 0;
  93     sc_.varstrs = nullptr;
  94     sc_.top     = nullptr;
  95     gmx_ana_index_clear(&sc_.gall);
  96     sc_.mempool = nullptr;
  97     sc_.symtab  = std::make_unique<SelectionParserSymbolTable>();
  98     gmx_ana_index_clear(&requiredAtoms_);
  99     gmx_ana_selmethod_register_defaults(sc_.symtab.get());
 100 }
 101
 102
 103 SelectionCollection::Impl::~Impl()
 104 {
 105     clearSymbolTable();
 106     // The tree must be freed before the SelectionData objects, since the
 107     // tree may hold references to the position data in SelectionData.
 108     sc_.root.reset();
 109     sc_.sel.clear();
 110     for (int i = 0; i < sc_.nvars; ++i)
 111     {
 112         sfree(sc_.varstrs[i]);
 113     }
 114     sfree(sc_.varstrs);
 115     gmx_ana_index_deinit(&sc_.gall);
 116     if (sc_.mempool)
 117     {
 118         _gmx_sel_mempool_destroy(sc_.mempool);
 119     }
 120     gmx_ana_index_deinit(&requiredAtoms_);
 121 }
 122
 123
 124 void SelectionCollection::Impl::clearSymbolTable()
 125 {
 126     sc_.symtab.reset();
 127 }
 128
 129
 130 namespace
 131 {
 132
 133 /*! \brief
 134  * Reads a single selection line from stdin.
 135  *
 136  * \param[in]  inputStream   Stream to read from (typically the StandardInputStream).
 137  * \param[in]  statusWriter  Stream to print prompts to (if NULL, no output is done).
 138  * \param[out] line          The read line in stored here.
 139  * \returns true if something was read, false if at end of input.
 140  *
 141  * Handles line continuation, reading also the continuing line(s) in one call.
 142  */
 143 bool promptLine(TextInputStream* inputStream, TextWriter* statusWriter, 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, 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     } while (status == YYPUSH_MORE);
 202     _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
 203     return status;
 204 }
 205
 206 /*! \brief
 207  * Print current status in response to empty line in interactive input.
 208  *
 209  * \param[in] writer         Writer to use for the output.
 210  * \param[in] sc             Selection collection data structure.
 211  * \param[in] grps           Available index groups.
 212  * \param[in] firstSelection Index of first selection from this interactive
 213  *     session.
 214  * \param[in] maxCount       Maximum number of selections.
 215  * \param[in] context        Context to print for what the selections are for.
 216  * \param[in] bFirst         Whether this is the header that is printed before
 217  *     any user input.
 218  *
 219  * Prints the available index groups and currently provided selections.
 220  */
 221 void printCurrentStatus(TextWriter*              writer,
 222                         gmx_ana_selcollection_t* sc,
 223                         gmx_ana_indexgrps_t*     grps,
 224                         size_t                   firstSelection,
 225                         int                      maxCount,
 226                         const std::string&       context,
 227                         bool                     bFirst)
 228 {
 229     if (grps != nullptr)
 230     {
 231         writer->writeLine("Available static index groups:");
 232         gmx_ana_indexgrps_print(writer, grps, 0);
 233     }
 234     writer->writeString("Specify ");
 235     if (maxCount < 0)
 236     {
 237         writer->writeString("any number of selections");
 238     }
 239     else if (maxCount == 1)
 240     {
 241         writer->writeString("a selection");
 242     }
 243     else
 244     {
 245         writer->writeString(formatString("%d selections", maxCount));
 246     }
 247     writer->writeString(formatString("%s%s:\n", context.empty() ? "" : " ", context.c_str()));
 248     writer->writeString(
 249             formatString("(one per line, <enter> for status/groups, 'help' for help%s)\n",
 250                          maxCount < 0 ? ", Ctrl-D to end" : ""));
 251     if (!bFirst && (sc->nvars > 0 || sc->sel.size() > firstSelection))
 252     {
 253         writer->writeLine("Currently provided selections:");
 254         for (int i = 0; i < sc->nvars; ++i)
 255         {
 256             writer->writeString(formatString("     %s\n", sc->varstrs[i]));
 257         }
 258         for (size_t i = firstSelection; i < sc->sel.size(); ++i)
 259         {
 260             writer->writeString(formatString(" %2d. %s\n", static_cast<int>(i - firstSelection + 1),
 261                                              sc->sel[i]->selectionText()));
 262         }
 263         if (maxCount > 0)
 264         {
 265             const int remaining = maxCount - static_cast<int>(sc->sel.size() - firstSelection);
 266             writer->writeString(formatString("(%d more selection%s required)\n", remaining,
 267                                              remaining > 1 ? "s" : ""));
 268         }
 269     }
 270 }
 271
 272 /*! \brief
 273  * Prints selection help in interactive selection input.
 274  *
 275  * \param[in] writer Writer to use for the output.
 276  * \param[in] sc    Selection collection data structure.
 277  * \param[in] line  Line of user input requesting help (starting with `help`).
 278  *
 279  * Initializes the selection help if not yet initialized, and finds the help
 280  * topic based on words on the input line.
 281  */
 282 void printHelp(TextWriter* writer, gmx_ana_selcollection_t* sc, const std::string& line)
 283 {
 284     if (sc->rootHelp.get() == nullptr)
 285     {
 286         sc->rootHelp = createSelectionHelpTopic();
 287     }
 288     HelpWriterContext context(writer, eHelpOutputFormat_Console);
 289     HelpManager       manager(*sc->rootHelp, context);
 290     try
 291     {
 292         std::vector<std::string>                 topic = splitString(line);
 293         std::vector<std::string>::const_iterator value;
 294         // First item in the list is the 'help' token.
 295         for (value = topic.begin() + 1; value != topic.end(); ++value)
 296         {
 297             manager.enterTopic(*value);
 298         }
 299     }
 300     catch (const InvalidInputError& ex)
 301     {
 302         writer->writeLine(ex.what());
 303         return;
 304     }
 305     manager.writeCurrentTopic();
 306 }
 307
 308 /*! \brief
 309  * Helper function that runs the parser once the tokenizer has been
 310  * initialized.
 311  *
 312  * \param[in,out] scanner       Scanner data structure.
 313  * \param[in]     inputStream   Stream to use for input (currently only with
 314  *      `bInteractive==true`).
 315  * \param[in]     bInteractive  Whether to use a line-based reading
 316  *      algorithm designed for interactive input.
 317  * \param[in]     maxnr   Maximum number of selections to parse
 318  *      (if -1, parse as many as provided by the user).
 319  * \param[in]     context Context to print for what the selections are for.
 320  * \returns       Vector of parsed selections.
 321  * \throws        std::bad_alloc if out of memory.
 322  * \throws        InvalidInputError if there is a parsing error.
 323  *
 324  * Used internally to implement parseInteractive(), parseFromFile() and
 325  * parseFromString().
 326  */
 327 SelectionList runParser(yyscan_t           scanner,
 328                         TextInputStream*   inputStream,
 329                         bool               bInteractive,
 330                         int                maxnr,
 331                         const std::string& context)
 332 {
 333     std::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         std::shared_ptr<_gmx_sel_yypstate> parserState(_gmx_sel_yypstate_new(), &_gmx_sel_yypstate_delete);
 340         if (bInteractive)
 341         {
 342             TextWriter* statusWriter = _gmx_sel_lexer_get_status_writer(scanner);
 343             if (statusWriter != nullptr)
 344             {
 345                 printCurrentStatus(statusWriter, sc, grps, oldCount, maxnr, context, true);
 346             }
 347             std::string line;
 348             int         status;
 349             while (promptLine(inputStream, statusWriter, &line))
 350             {
 351                 if (statusWriter != nullptr)
 352                 {
 353                     line = stripString(line);
 354                     if (line.empty())
 355                     {
 356                         printCurrentStatus(statusWriter, sc, grps, oldCount, maxnr, context, false);
 357                         continue;
 358                     }
 359                     if (startsWith(line, "help") && (line[4] == 0 || (std::isspace(line[4]) != 0)))
 360                     {
 361                         printHelp(statusWriter, sc, line);
 362                         continue;
 363                     }
 364                 }
 365                 line.append("\n");
 366                 _gmx_sel_set_lex_input_str(scanner, line.c_str());
 367                 status = runParserLoop(scanner, parserState.get(), true);
 368                 if (status != YYPUSH_MORE)
 369                 {
 370                     // TODO: Check if there is more input, and issue an
 371                     // error/warning if some input was ignored.
 372                     goto early_termination;
 373                 }
 374             }
 375             {
 376                 YYLTYPE location;
 377                 status = _gmx_sel_yypush_parse(parserState.get(), 0, nullptr, &location, scanner);
 378             }
 379             // TODO: Remove added selections from the collection if parsing failed?
 380             _gmx_sel_lexer_rethrow_exception_if_occurred(scanner);
 381         early_termination:
 382             GMX_RELEASE_ASSERT(status == 0, "Parser errors should have resulted in an exception");
 383         }
 384         else
 385         {
 386             int status = runParserLoop(scanner, parserState.get(), false);
 387             GMX_RELEASE_ASSERT(status == 0, "Parser errors should have resulted in an exception");
 388         }
 389     }
 390     scannerGuard.reset();
 391     int nr = sc->sel.size() - oldCount;
 392     if (maxnr > 0 && nr != maxnr)
 393     {
 394         std::string message = formatString("Too few selections provided; got %d, expected %d", nr, maxnr);
 395         GMX_THROW(InvalidInputError(message));
 396     }
 397
 398     SelectionList                     result;
 399     SelectionDataList::const_iterator i;
 400     result.reserve(nr);
 401     for (i = sc->sel.begin() + oldCount; i != sc->sel.end(); ++i)
 402     {
 403         result.emplace_back(i->get());
 404     }
 405     return result;
 406 }
 407
 408 /*! \brief
 409  * Checks that index groups have valid atom indices.
 410  *
 411  * \param[in]    root    Root of selection tree to process.
 412  * \param[in]    natoms  Maximum number of atoms that the selections are set
 413  *     to evaluate.
 414  * \param        errors  Object for reporting any error messages.
 415  * \throws std::bad_alloc if out of memory.
 416  *
 417  * Recursively checks the selection tree for index groups.
 418  * Each found group is checked that it only contains atom indices that match
 419  * the topology/maximum number of atoms set for the selection collection.
 420  * Any issues are reported to \p errors.
 421  */
 422 void checkExternalGroups(const SelectionTreeElementPointer& root, int natoms, ExceptionInitializer* errors)
 423 {
 424     if (root->type == SEL_CONST && root->v.type == GROUP_VALUE)
 425     {
 426         try
 427         {
 428             root->checkIndexGroup(natoms);
 429         }
 430         catch (const UserInputError&)
 431         {
 432             errors->addCurrentExceptionAsNested();
 433         }
 434     }
 435
 436     SelectionTreeElementPointer child = root->child;
 437     while (child)
 438     {
 439         checkExternalGroups(child, natoms, errors);
 440         child = child->next;
 441     }
 442 }
 443
 444 //! Checks whether the given topology properties are available.
 445 void checkTopologyProperties(const gmx_mtop_t* top, const SelectionTopologyProperties& props)
 446 {
 447     if (top == nullptr)
 448     {
 449         if (props.hasAny())
 450         {
 451             GMX_THROW(InconsistentInputError(
 452                     "Selection requires topology information, but none provided"));
 453         }
 454         return;
 455     }
 456     if (props.needsMasses && !gmx_mtop_has_masses(top))
 457     {
 458         GMX_THROW(InconsistentInputError(
 459                 "Selection requires mass information, but it is not available in the topology"));
 460     }
 461 }
 462
 463 } // namespace
 464
 465
 466 void SelectionCollection::Impl::resolveExternalGroups(const SelectionTreeElementPointer& root,
 467                                                       ExceptionInitializer*              errors)
 468 {
 469
 470     if (root->type == SEL_GROUPREF)
 471     {
 472         try
 473         {
 474             root->resolveIndexGroupReference(grps_, sc_.gall.isize);
 475         }
 476         catch (const UserInputError&)
 477         {
 478             errors->addCurrentExceptionAsNested();
 479         }
 480     }
 481
 482     SelectionTreeElementPointer child = root->child;
 483     while (child)
 484     {
 485         resolveExternalGroups(child, errors);
 486         root->flags |= (child->flags & SEL_UNSORTED);
 487         child = child->next;
 488     }
 489 }
 490
 491
 492 bool SelectionCollection::Impl::areForcesRequested() const
 493 {
 494     for (const auto& sel : sc_.sel)
 495     {
 496         if (sel->hasFlag(gmx::efSelection_EvaluateForces))
 497         {
 498             return true;
 499         }
 500     }
 501     return false;
 502 }
 503
 504
 505 SelectionTopologyProperties
 506 SelectionCollection::Impl::requiredTopologyPropertiesForPositionType(const std::string& post, bool forces) const
 507 {
 508     SelectionTopologyProperties props;
 509     if (!post.empty())
 510     {
 511         switch (PositionCalculationCollection::requiredTopologyInfoForType(post.c_str(), forces))
 512         {
 513             case PositionCalculationCollection::RequiredTopologyInfo::None: break;
 514             case PositionCalculationCollection::RequiredTopologyInfo::Topology:
 515                 props.merge(SelectionTopologyProperties::topology());
 516                 break;
 517             case PositionCalculationCollection::RequiredTopologyInfo::TopologyAndMasses:
 518                 props.merge(SelectionTopologyProperties::masses());
 519                 break;
 520         }
 521     }
 522     return props;
 523 }
 524
 525
 526 /********************************************************************
 527  * SelectionCollection
 528  */
 529
 530 SelectionCollection::SelectionCollection() : impl_(new Impl) {}
 531
 532
 533 SelectionCollection::~SelectionCollection() {}
 534
 535
 536 void SelectionCollection::initOptions(IOptionsContainer* options, SelectionTypeOption selectionTypeOption)
 537 {
 538     static const EnumerationArray<Impl::DebugLevel, const char*> s_debugLevelNames = {
 539         { "no", "basic", "compile", "eval", "full" }
 540     };
 541
 542     const char* const* postypes = PositionCalculationCollection::typeEnumValues;
 543     options->addOption(StringOption("selrpos")
 544                                .enumValueFromNullTerminatedArray(postypes)
 545                                .store(&impl_->rpost_)
 546                                .defaultValue(postypes[0])
 547                                .description("Selection reference positions"));
 548     if (selectionTypeOption == IncludeSelectionTypeOption)
 549     {
 550         options->addOption(StringOption("seltype")
 551                                    .enumValueFromNullTerminatedArray(postypes)
 552                                    .store(&impl_->spost_)
 553                                    .defaultValue(postypes[0])
 554                                    .description("Default selection output positions"));
 555     }
 556     else
 557     {
 558         impl_->spost_ = postypes[0];
 559     }
 560     GMX_RELEASE_ASSERT(impl_->debugLevel_ != Impl::DebugLevel::Count, "Debug level out of range");
 561     options->addOption(EnumOption<Impl::DebugLevel>("seldebug")
 562                                .hidden(impl_->debugLevel_ == Impl::DebugLevel::None)
 563                                .enumValue(s_debugLevelNames)
 564                                .store(&impl_->debugLevel_)
 565                                .description("Print out selection trees for debugging"));
 566 }
 567
 568
 569 void SelectionCollection::setReferencePosType(const char* type)
 570 {
 571     GMX_RELEASE_ASSERT(type != nullptr, "Cannot assign NULL position type");
 572     // Check that the type is valid, throw if it is not.
 573     e_poscalc_t dummytype;
 574     int         dummyflags;
 575     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
 576     impl_->rpost_ = type;
 577 }
 578
 579
 580 void SelectionCollection::setOutputPosType(const char* type)
 581 {
 582     GMX_RELEASE_ASSERT(type != nullptr, "Cannot assign NULL position type");
 583     // Check that the type is valid, throw if it is not.
 584     e_poscalc_t dummytype;
 585     int         dummyflags;
 586     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
 587     impl_->spost_ = type;
 588 }
 589
 590
 591 void SelectionCollection::setDebugLevel(int debugLevel)
 592 {
 593     impl_->debugLevel_ = Impl::DebugLevel(debugLevel);
 594 }
 595
 596
 597 void SelectionCollection::setTopology(gmx_mtop_t* top, int natoms)
 598 {
 599     GMX_RELEASE_ASSERT(natoms > 0 || top != nullptr,
 600                        "The number of atoms must be given if there is no topology");
 601     checkTopologyProperties(top, requiredTopologyProperties());
 602     // Get the number of atoms from the topology if it is not given.
 603     if (natoms <= 0)
 604     {
 605         natoms = top->natoms;
 606     }
 607     if (impl_->bExternalGroupsSet_)
 608     {
 609         ExceptionInitializer        errors("Invalid index group references encountered");
 610         SelectionTreeElementPointer root = impl_->sc_.root;
 611         while (root)
 612         {
 613             checkExternalGroups(root, natoms, &errors);
 614             root = root->next;
 615         }
 616         if (errors.hasNestedExceptions())
 617         {
 618             GMX_THROW(InconsistentInputError(errors));
 619         }
 620     }
 621     gmx_ana_selcollection_t* sc = &impl_->sc_;
 622     // Do this first, as it allocates memory, while the others don't throw.
 623     gmx_ana_index_init_simple(&sc->gall, natoms);
 624     sc->top = top;
 625     sc->pcc.setTopology(top);
 626 }
 627
 628
 629 void SelectionCollection::setIndexGroups(gmx_ana_indexgrps_t* grps)
 630 {
 631     GMX_RELEASE_ASSERT(grps == nullptr || !impl_->bExternalGroupsSet_,
 632                        "Can only set external groups once or clear them afterwards");
 633     impl_->grps_               = grps;
 634     impl_->bExternalGroupsSet_ = true;
 635
 636     ExceptionInitializer        errors("Invalid index group reference(s)");
 637     SelectionTreeElementPointer root = impl_->sc_.root;
 638     while (root)
 639     {
 640         impl_->resolveExternalGroups(root, &errors);
 641         root->checkUnsortedAtoms(true, &errors);
 642         root = root->next;
 643     }
 644     if (errors.hasNestedExceptions())
 645     {
 646         GMX_THROW(InconsistentInputError(errors));
 647     }
 648     for (size_t i = 0; i < impl_->sc_.sel.size(); ++i)
 649     {
 650         impl_->sc_.sel[i]->refreshName();
 651     }
 652 }
 653
 654 SelectionTopologyProperties SelectionCollection::requiredTopologyProperties() const
 655 {
 656     SelectionTopologyProperties props;
 657
 658     // These should not throw, because has been checked earlier.
 659     props.merge(impl_->requiredTopologyPropertiesForPositionType(impl_->rpost_, false));
 660     const bool forcesRequested = impl_->areForcesRequested();
 661     props.merge(impl_->requiredTopologyPropertiesForPositionType(impl_->spost_, forcesRequested));
 662
 663     SelectionTreeElementPointer sel = impl_->sc_.root;
 664     while (sel && !props.hasAll())
 665     {
 666         props.merge(sel->requiredTopologyProperties());
 667         sel = sel->next;
 668     }
 669     return props;
 670 }
 671
 672
 673 bool SelectionCollection::requiresIndexGroups() const
 674 {
 675     SelectionTreeElementPointer sel = impl_->sc_.root;
 676     while (sel)
 677     {
 678         if (sel->requiresIndexGroups())
 679         {
 680             return true;
 681         }
 682         sel = sel->next;
 683     }
 684     return false;
 685 }
 686
 687
 688 SelectionList SelectionCollection::parseFromStdin(int count, bool bInteractive, const std::string& context)
 689 {
 690     return parseInteractive(count, &StandardInputStream::instance(),
 691                             bInteractive ? &TextOutputFile::standardError() : nullptr, context);
 692 }
 693
 694 namespace
 695 {
 696
 697 //! Helper function to initialize status writer for interactive selection parsing.
 698 std::unique_ptr<TextWriter> initStatusWriter(TextOutputStream* statusStream)
 699 {
 700     std::unique_ptr<TextWriter> statusWriter;
 701     if (statusStream != nullptr)
 702     {
 703         statusWriter = std::make_unique<TextWriter>(statusStream);
 704         statusWriter->wrapperSettings().setLineLength(78);
 705     }
 706     return statusWriter;
 707 }
 708
 709 } // namespace
 710
 711 SelectionList SelectionCollection::parseInteractive(int                count,
 712                                                     TextInputStream*   inputStream,
 713                                                     TextOutputStream*  statusStream,
 714                                                     const std::string& context)
 715 {
 716     yyscan_t scanner;
 717
 718     const std::unique_ptr<TextWriter> statusWriter(initStatusWriter(statusStream));
 719     _gmx_sel_init_lexer(&scanner, &impl_->sc_, statusWriter.get(), count,
 720                         impl_->bExternalGroupsSet_, impl_->grps_);
 721     return runParser(scanner, inputStream, true, count, context);
 722 }
 723
 724
 725 SelectionList SelectionCollection::parseFromFile(const std::string& filename)
 726 {
 727
 728     try
 729     {
 730         yyscan_t      scanner;
 731         TextInputFile file(filename);
 732         // TODO: Exception-safe way of using the lexer.
 733         _gmx_sel_init_lexer(&scanner, &impl_->sc_, nullptr, -1, impl_->bExternalGroupsSet_, impl_->grps_);
 734         _gmx_sel_set_lex_input_file(scanner, file.handle());
 735         return runParser(scanner, nullptr, false, -1, std::string());
 736     }
 737     catch (GromacsException& ex)
 738     {
 739         ex.prependContext(formatString("Error in parsing selections from file '%s'", filename.c_str()));
 740         throw;
 741     }
 742 }
 743
 744
 745 SelectionList SelectionCollection::parseFromString(const std::string& str)
 746 {
 747     yyscan_t scanner;
 748
 749     _gmx_sel_init_lexer(&scanner, &impl_->sc_, nullptr, -1, impl_->bExternalGroupsSet_, impl_->grps_);
 750     _gmx_sel_set_lex_input_str(scanner, str.c_str());
 751     return runParser(scanner, nullptr, false, -1, std::string());
 752 }
 753
 754
 755 void SelectionCollection::compile()
 756 {
 757     checkTopologyProperties(impl_->sc_.top, requiredTopologyProperties());
 758     if (!impl_->bExternalGroupsSet_)
 759     {
 760         setIndexGroups(nullptr);
 761     }
 762     if (impl_->debugLevel_ != Impl::DebugLevel::None)
 763     {
 764         printTree(stderr, false);
 765     }
 766
 767     SelectionCompiler compiler;
 768     compiler.compile(this);
 769
 770     if (impl_->debugLevel_ != Impl::DebugLevel::None)
 771     {
 772         std::fprintf(stderr, "\n");
 773         printTree(stderr, false);
 774         std::fprintf(stderr, "\n");
 775         impl_->sc_.pcc.printTree(stderr);
 776         std::fprintf(stderr, "\n");
 777     }
 778     impl_->sc_.pcc.initEvaluation();
 779     if (impl_->debugLevel_ != Impl::DebugLevel::None)
 780     {
 781         impl_->sc_.pcc.printTree(stderr);
 782         std::fprintf(stderr, "\n");
 783     }
 784
 785     // TODO: It would be nicer to associate the name of the selection option
 786     // (if available) to the error message.
 787     SelectionDataList::const_iterator iter;
 788     for (iter = impl_->sc_.sel.begin(); iter != impl_->sc_.sel.end(); ++iter)
 789     {
 790         const internal::SelectionData& sel = **iter;
 791         if (sel.hasFlag(efSelection_OnlyAtoms))
 792         {
 793             if (!sel.hasOnlyAtoms())
 794             {
 795                 std::string message = formatString(
 796                         "Selection '%s' does not evaluate to individual atoms. "
 797                         "This is not allowed in this context.",
 798                         sel.selectionText());
 799                 GMX_THROW(InvalidInputError(message));
 800             }
 801             if (sel.hasFlag(efSelection_OnlySorted))
 802             {
 803                 if (!sel.hasSortedAtomIndices())
 804                 {
 805                     const std::string message = formatString(
 806                             "Selection '%s' does not evaluate to atoms in an "
 807                             "ascending (sorted) order. "
 808                             "This is not allowed in this context.",
 809                             sel.selectionText());
 810                     GMX_THROW(InvalidInputError(message));
 811                 }
 812             }
 813         }
 814         if (sel.hasFlag(efSelection_DisallowEmpty))
 815         {
 816             if (sel.posCount() == 0)
 817             {
 818                 std::string message =
 819                         formatString("Selection '%s' never matches any atoms.", sel.selectionText());
 820                 GMX_THROW(InvalidInputError(message));
 821             }
 822         }
 823     }
 824     impl_->rpost_.clear();
 825     impl_->spost_.clear();
 826 }
 827
 828
 829 void SelectionCollection::evaluate(t_trxframe* fr, t_pbc* pbc)
 830 {
 831     checkTopologyProperties(impl_->sc_.top, requiredTopologyProperties());
 832     if (fr->bIndex)
 833     {
 834         gmx_ana_index_t g;
 835         gmx_ana_index_set(&g, fr->natoms, fr->index, 0);
 836         GMX_RELEASE_ASSERT(gmx_ana_index_check_sorted(&g),
 837                            "Only trajectories with atoms in ascending order "
 838                            "are currently supported");
 839         if (!gmx_ana_index_contains(&g, &impl_->requiredAtoms_))
 840         {
 841             const std::string message = formatString(
 842                     "Trajectory does not contain all atoms required for "
 843                     "evaluating the provided selections.");
 844             GMX_THROW(InconsistentInputError(message));
 845         }
 846     }
 847     else
 848     {
 849         const int maxAtomIndex = gmx_ana_index_get_max_index(&impl_->requiredAtoms_);
 850         if (fr->natoms <= maxAtomIndex)
 851         {
 852             const std::string message = formatString(
 853                     "Trajectory has less atoms (%d) than what is required for "
 854                     "evaluating the provided selections (atoms up to index %d "
 855                     "are required).",
 856                     fr->natoms, maxAtomIndex + 1);
 857             GMX_THROW(InconsistentInputError(message));
 858         }
 859     }
 860     impl_->sc_.pcc.initFrame(fr);
 861
 862     SelectionEvaluator evaluator;
 863     evaluator.evaluate(this, fr, pbc);
 864
 865     if (impl_->debugLevel_ == Impl::DebugLevel::Evaluated || impl_->debugLevel_ == Impl::DebugLevel::Full)
 866     {
 867         std::fprintf(stderr, "\n");
 868         printTree(stderr, true);
 869     }
 870 }
 871
 872
 873 void SelectionCollection::evaluateFinal(int nframes)
 874 {
 875     SelectionEvaluator evaluator;
 876     evaluator.evaluateFinal(this, nframes);
 877 }
 878
 879
 880 void SelectionCollection::printTree(FILE* fp, bool bValues) const
 881 {
 882     SelectionTreeElementPointer sel = impl_->sc_.root;
 883     while (sel)
 884     {
 885         _gmx_selelem_print_tree(fp, *sel, bValues, 0);
 886         sel = sel->next;
 887     }
 888 }
 889
 890
 891 void SelectionCollection::printXvgrInfo(FILE* out) const
 892 {
 893     const gmx_ana_selcollection_t& sc = impl_->sc_;
 894     std::fprintf(out, "# Selections:\n");
 895     for (int i = 0; i < sc.nvars; ++i)
 896     {
 897         std::fprintf(out, "#   %s\n", sc.varstrs[i]);
 898     }
 899     for (size_t i = 0; i < sc.sel.size(); ++i)
 900     {
 901         std::fprintf(out, "#   %s\n", sc.sel[i]->selectionText());
 902     }
 903     std::fprintf(out, "#\n");
 904 }
 905
 906 } // namespace gmx