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