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),
  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(
 261                     " %2d. %s\n", static_cast<int>(i - firstSelection + 1), 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(
 267                     "(%d more selection%s required)\n", remaining, 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)
 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() = default;
 534
 535 SelectionCollection::SelectionCollection(const SelectionCollection& rhs) :
 536     impl_(std::make_unique<Impl>())
 537 {
 538     setReferencePosType(rhs.impl_->rpost_.empty() ? PositionCalculationCollection::typeEnumValues[0]
 539                                                   : rhs.impl_->rpost_.c_str());
 540     setOutputPosType(rhs.impl_->spost_.empty() ? PositionCalculationCollection::typeEnumValues[0]
 541                                                : rhs.impl_->spost_.c_str());
 542     setDebugLevel(static_cast<int>(rhs.impl_->debugLevel_));
 543
 544     for (size_t i = 0; i < rhs.impl_->sc_.sel.size(); i++)
 545     {
 546         const auto& selectionOption = rhs.impl_->sc_.sel[i];
 547         parseFromString(selectionOption->selectionText());
 548         impl_->sc_.sel[i]->setFlags(selectionOption->flags());
 549     }
 550
 551     // Topology has been initialized in rhs if top is non-null or natoms is set.
 552     // Note this needs to be set after selections are parsed to register topology requirements properly.
 553     if (rhs.impl_->sc_.top != nullptr || rhs.impl_->sc_.gall.isize > 0)
 554     {
 555         setTopology(rhs.impl_->sc_.top, rhs.impl_->sc_.gall.isize);
 556         gmx_ana_index_copy(&impl_->sc_.gall, &rhs.impl_->sc_.gall, /*balloc=*/false);
 557     }
 558
 559     if (rhs.impl_->grps_ != nullptr)
 560     {
 561         setIndexGroups(rhs.impl_->grps_);
 562     }
 563
 564     // Only compile the selection if rhs is compiled.
 565     if (rhs.impl_->sc_.mempool != nullptr)
 566     {
 567         compile();
 568     }
 569 }
 570
 571 SelectionCollection& SelectionCollection::operator=(SelectionCollection rhs)
 572 {
 573     rhs.swap(*this);
 574     return *this;
 575 }
 576
 577 void SelectionCollection::swap(SelectionCollection& rhs)
 578 {
 579     using std::swap;
 580     swap(impl_, rhs.impl_);
 581 }
 582
 583 void SelectionCollection::initOptions(IOptionsContainer* options, SelectionTypeOption selectionTypeOption)
 584 {
 585     static const EnumerationArray<Impl::DebugLevel, const char*> s_debugLevelNames = {
 586         { "no", "basic", "compile", "eval", "full" }
 587     };
 588
 589     const char* const* postypes = PositionCalculationCollection::typeEnumValues;
 590     options->addOption(StringOption("selrpos")
 591                                .enumValueFromNullTerminatedArray(postypes)
 592                                .store(&impl_->rpost_)
 593                                .defaultValue(postypes[0])
 594                                .description("Selection reference positions"));
 595     if (selectionTypeOption == IncludeSelectionTypeOption)
 596     {
 597         options->addOption(StringOption("seltype")
 598                                    .enumValueFromNullTerminatedArray(postypes)
 599                                    .store(&impl_->spost_)
 600                                    .defaultValue(postypes[0])
 601                                    .description("Default selection output positions"));
 602     }
 603     else
 604     {
 605         impl_->spost_ = postypes[0];
 606     }
 607     GMX_RELEASE_ASSERT(impl_->debugLevel_ != Impl::DebugLevel::Count, "Debug level out of range");
 608     options->addOption(EnumOption<Impl::DebugLevel>("seldebug")
 609                                .hidden(impl_->debugLevel_ == Impl::DebugLevel::None)
 610                                .enumValue(s_debugLevelNames)
 611                                .store(&impl_->debugLevel_)
 612                                .description("Print out selection trees for debugging"));
 613 }
 614
 615
 616 void SelectionCollection::setReferencePosType(const char* type)
 617 {
 618     GMX_RELEASE_ASSERT(type != nullptr, "Cannot assign NULL position type");
 619     // Check that the type is valid, throw if it is not.
 620     e_poscalc_t dummytype;
 621     int         dummyflags;
 622     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
 623     impl_->rpost_ = type;
 624 }
 625
 626
 627 void SelectionCollection::setOutputPosType(const char* type)
 628 {
 629     GMX_RELEASE_ASSERT(type != nullptr, "Cannot assign NULL position type");
 630     // Check that the type is valid, throw if it is not.
 631     e_poscalc_t dummytype;
 632     int         dummyflags;
 633     PositionCalculationCollection::typeFromEnum(type, &dummytype, &dummyflags);
 634     impl_->spost_ = type;
 635 }
 636
 637
 638 void SelectionCollection::setDebugLevel(int debugLevel)
 639 {
 640     impl_->debugLevel_ = Impl::DebugLevel(debugLevel);
 641 }
 642
 643
 644 void SelectionCollection::setTopology(const gmx_mtop_t* top, int natoms)
 645 {
 646     GMX_RELEASE_ASSERT(natoms > 0 || top != nullptr,
 647                        "The number of atoms must be given if there is no topology");
 648     checkTopologyProperties(top, requiredTopologyProperties());
 649     // Get the number of atoms from the topology if it is not given.
 650     if (natoms <= 0)
 651     {
 652         natoms = top->natoms;
 653     }
 654     if (impl_->bExternalGroupsSet_)
 655     {
 656         ExceptionInitializer        errors("Invalid index group references encountered");
 657         SelectionTreeElementPointer root = impl_->sc_.root;
 658         while (root)
 659         {
 660             checkExternalGroups(root, natoms, &errors);
 661             root = root->next;
 662         }
 663         if (errors.hasNestedExceptions())
 664         {
 665             GMX_THROW(InconsistentInputError(errors));
 666         }
 667     }
 668     gmx_ana_selcollection_t* sc = &impl_->sc_;
 669     // Do this first, as it allocates memory, while the others don't throw.
 670     gmx_ana_index_init_simple(&sc->gall, natoms);
 671     sc->top = top;
 672     sc->pcc.setTopology(top);
 673 }
 674
 675
 676 void SelectionCollection::setIndexGroups(gmx_ana_indexgrps_t* grps)
 677 {
 678     GMX_RELEASE_ASSERT(grps == nullptr || !impl_->bExternalGroupsSet_,
 679                        "Can only set external groups once or clear them afterwards");
 680     impl_->grps_               = grps;
 681     impl_->bExternalGroupsSet_ = true;
 682
 683     ExceptionInitializer        errors("Invalid index group reference(s)");
 684     SelectionTreeElementPointer root = impl_->sc_.root;
 685     while (root)
 686     {
 687         impl_->resolveExternalGroups(root, &errors);
 688         root->checkUnsortedAtoms(true, &errors);
 689         root = root->next;
 690     }
 691     if (errors.hasNestedExceptions())
 692     {
 693         GMX_THROW(InconsistentInputError(errors));
 694     }
 695     for (size_t i = 0; i < impl_->sc_.sel.size(); ++i)
 696     {
 697         impl_->sc_.sel[i]->refreshName();
 698     }
 699 }
 700
 701 SelectionTopologyProperties SelectionCollection::requiredTopologyProperties() const
 702 {
 703     SelectionTopologyProperties props;
 704
 705     // These should not throw, because has been checked earlier.
 706     props.merge(impl_->requiredTopologyPropertiesForPositionType(impl_->rpost_, false));
 707     const bool forcesRequested = impl_->areForcesRequested();
 708     props.merge(impl_->requiredTopologyPropertiesForPositionType(impl_->spost_, forcesRequested));
 709
 710     SelectionTreeElementPointer sel = impl_->sc_.root;
 711     while (sel && !props.hasAll())
 712     {
 713         props.merge(sel->requiredTopologyProperties());
 714         sel = sel->next;
 715     }
 716     return props;
 717 }
 718
 719
 720 bool SelectionCollection::requiresIndexGroups() const
 721 {
 722     SelectionTreeElementPointer sel = impl_->sc_.root;
 723     while (sel)
 724     {
 725         if (sel->requiresIndexGroups())
 726         {
 727             return true;
 728         }
 729         sel = sel->next;
 730     }
 731     return false;
 732 }
 733
 734
 735 SelectionList SelectionCollection::parseFromStdin(int count, bool bInteractive, const std::string& context)
 736 {
 737     StandardInputStream inputStream;
 738     return parseInteractive(
 739             count, &inputStream, bInteractive ? &TextOutputFile::standardError() : nullptr, context);
 740 }
 741
 742 namespace
 743 {
 744
 745 //! Helper function to initialize status writer for interactive selection parsing.
 746 std::unique_ptr<TextWriter> initStatusWriter(TextOutputStream* statusStream)
 747 {
 748     std::unique_ptr<TextWriter> statusWriter;
 749     if (statusStream != nullptr)
 750     {
 751         statusWriter = std::make_unique<TextWriter>(statusStream);
 752         statusWriter->wrapperSettings().setLineLength(78);
 753     }
 754     return statusWriter;
 755 }
 756
 757 } // namespace
 758
 759 SelectionList SelectionCollection::parseInteractive(int                count,
 760                                                     TextInputStream*   inputStream,
 761                                                     TextOutputStream*  statusStream,
 762                                                     const std::string& context)
 763 {
 764     yyscan_t scanner;
 765
 766     const std::unique_ptr<TextWriter> statusWriter(initStatusWriter(statusStream));
 767     _gmx_sel_init_lexer(
 768             &scanner, &impl_->sc_, statusWriter.get(), count, impl_->bExternalGroupsSet_, impl_->grps_);
 769     return runParser(scanner, inputStream, true, count, context);
 770 }
 771
 772
 773 SelectionList SelectionCollection::parseFromFile(const std::string& filename)
 774 {
 775
 776     try
 777     {
 778         yyscan_t      scanner;
 779         TextInputFile file(filename);
 780         // TODO: Exception-safe way of using the lexer.
 781         _gmx_sel_init_lexer(&scanner, &impl_->sc_, nullptr, -1, impl_->bExternalGroupsSet_, impl_->grps_);
 782         _gmx_sel_set_lex_input_file(scanner, file.handle());
 783         return runParser(scanner, nullptr, false, -1, std::string());
 784     }
 785     catch (GromacsException& ex)
 786     {
 787         ex.prependContext(formatString("Error in parsing selections from file '%s'", filename.c_str()));
 788         throw;
 789     }
 790 }
 791
 792
 793 SelectionList SelectionCollection::parseFromString(const std::string& str)
 794 {
 795     yyscan_t scanner;
 796
 797     _gmx_sel_init_lexer(&scanner, &impl_->sc_, nullptr, -1, impl_->bExternalGroupsSet_, impl_->grps_);
 798     _gmx_sel_set_lex_input_str(scanner, str.c_str());
 799     return runParser(scanner, nullptr, false, -1, std::string());
 800 }
 801
 802
 803 void SelectionCollection::compile()
 804 {
 805     checkTopologyProperties(impl_->sc_.top, requiredTopologyProperties());
 806     if (!impl_->bExternalGroupsSet_)
 807     {
 808         setIndexGroups(nullptr);
 809     }
 810     if (impl_->debugLevel_ != Impl::DebugLevel::None)
 811     {
 812         printTree(stderr, false);
 813     }
 814
 815     compileSelection(this);
 816
 817     if (impl_->debugLevel_ != Impl::DebugLevel::None)
 818     {
 819         std::fprintf(stderr, "\n");
 820         printTree(stderr, false);
 821         std::fprintf(stderr, "\n");
 822         impl_->sc_.pcc.printTree(stderr);
 823         std::fprintf(stderr, "\n");
 824     }
 825     impl_->sc_.pcc.initEvaluation();
 826     if (impl_->debugLevel_ != Impl::DebugLevel::None)
 827     {
 828         impl_->sc_.pcc.printTree(stderr);
 829         std::fprintf(stderr, "\n");
 830     }
 831
 832     // TODO: It would be nicer to associate the name of the selection option
 833     // (if available) to the error message.
 834     SelectionDataList::const_iterator iter;
 835     for (iter = impl_->sc_.sel.begin(); iter != impl_->sc_.sel.end(); ++iter)
 836     {
 837         const internal::SelectionData& sel = **iter;
 838         if (sel.hasFlag(efSelection_OnlyAtoms))
 839         {
 840             if (!sel.hasOnlyAtoms())
 841             {
 842                 std::string message = formatString(
 843                         "Selection '%s' does not evaluate to individual atoms. "
 844                         "This is not allowed in this context.",
 845                         sel.selectionText());
 846                 GMX_THROW(InvalidInputError(message));
 847             }
 848             if (sel.hasFlag(efSelection_OnlySorted))
 849             {
 850                 if (!sel.hasSortedAtomIndices())
 851                 {
 852                     const std::string message = formatString(
 853                             "Selection '%s' does not evaluate to atoms in an "
 854                             "ascending (sorted) order. "
 855                             "This is not allowed in this context.",
 856                             sel.selectionText());
 857                     GMX_THROW(InvalidInputError(message));
 858                 }
 859             }
 860         }
 861         if (sel.hasFlag(efSelection_DisallowEmpty))
 862         {
 863             if (sel.posCount() == 0)
 864             {
 865                 std::string message =
 866                         formatString("Selection '%s' never matches any atoms.", sel.selectionText());
 867                 GMX_THROW(InvalidInputError(message));
 868             }
 869         }
 870     }
 871     impl_->rpost_.clear();
 872     impl_->spost_.clear();
 873 }
 874
 875
 876 void SelectionCollection::evaluate(t_trxframe* fr, t_pbc* pbc)
 877 {
 878     checkTopologyProperties(impl_->sc_.top, requiredTopologyProperties());
 879     if (fr->bIndex)
 880     {
 881         gmx_ana_index_t g;
 882         gmx_ana_index_set(&g, fr->natoms, fr->index, 0);
 883         GMX_RELEASE_ASSERT(gmx_ana_index_check_sorted(&g),
 884                            "Only trajectories with atoms in ascending order "
 885                            "are currently supported");
 886         if (!gmx_ana_index_contains(&g, &impl_->requiredAtoms_))
 887         {
 888             const std::string message = formatString(
 889                     "Trajectory does not contain all atoms required for "
 890                     "evaluating the provided selections.");
 891             GMX_THROW(InconsistentInputError(message));
 892         }
 893     }
 894     else
 895     {
 896         const int maxAtomIndex = gmx_ana_index_get_max_index(&impl_->requiredAtoms_);
 897         if (fr->natoms <= maxAtomIndex)
 898         {
 899             const std::string message = formatString(
 900                     "Trajectory has less atoms (%d) than what is required for "
 901                     "evaluating the provided selections (atoms up to index %d "
 902                     "are required).",
 903                     fr->natoms,
 904                     maxAtomIndex + 1);
 905             GMX_THROW(InconsistentInputError(message));
 906         }
 907     }
 908     impl_->sc_.pcc.initFrame(fr);
 909
 910     SelectionEvaluator evaluator;
 911     evaluator.evaluate(this, fr, pbc);
 912
 913     if (impl_->debugLevel_ == Impl::DebugLevel::Evaluated || impl_->debugLevel_ == Impl::DebugLevel::Full)
 914     {
 915         std::fprintf(stderr, "\n");
 916         printTree(stderr, true);
 917     }
 918 }
 919
 920
 921 void SelectionCollection::evaluateFinal(int nframes)
 922 {
 923     SelectionEvaluator evaluator;
 924     evaluator.evaluateFinal(this, nframes);
 925 }
 926
 927 std::optional<Selection> SelectionCollection::selection(std::string_view selName) const
 928 {
 929     const auto& selections = impl_->sc_.sel;
 930     if (const auto foundIter = std::find_if(
 931                 selections.cbegin(),
 932                 selections.cend(),
 933                 [selName](const auto& selection) { return selection->name() == selName; });
 934         foundIter != selections.end())
 935     {
 936         return Selection(foundIter->get());
 937     }
 938     return std::nullopt;
 939 }
 940
 941
 942 void SelectionCollection::printTree(FILE* fp, bool bValues) const
 943 {
 944     SelectionTreeElementPointer sel = impl_->sc_.root;
 945     while (sel)
 946     {
 947         _gmx_selelem_print_tree(fp, *sel, bValues, 0);
 948         sel = sel->next;
 949     }
 950 }
 951
 952
 953 void SelectionCollection::printXvgrInfo(FILE* out) const
 954 {
 955     const gmx_ana_selcollection_t& sc = impl_->sc_;
 956     std::fprintf(out, "# Selections:\n");
 957     for (int i = 0; i < sc.nvars; ++i)
 958     {
 959         std::fprintf(out, "#   %s\n", sc.varstrs[i]);
 960     }
 961     for (size_t i = 0; i < sc.sel.size(); ++i)
 962     {
 963         std::fprintf(out, "#   %s\n", sc.sel[i]->selectionText());
 964     }
 965     std::fprintf(out, "#\n");
 966 }
 967
 968 void swap(SelectionCollection& lhs, SelectionCollection& rhs)
 969 {
 970     lhs.swap(rhs);
 971 }
 972
 973 } // namespace gmx