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