src/gromacs/analysisdata/modules/plot.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 classes in plot.h.
  39  *
  40  * \ingroup module_analysisdata
  41  * \author Teemu Murtola <teemu.murtola@gmail.com>
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "plot.h"
  46
  47 #include <cstdio>
  48 #include <cstring>
  49
  50 #include <string>
  51 #include <vector>
  52
  53 #include "gromacs/analysisdata/dataframe.h"
  54 #include "gromacs/fileio/gmxfio.h"
  55 #include "gromacs/fileio/oenv.h"
  56 #include "gromacs/fileio/xvgr.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/options/basicoptions.h"
  59 #include "gromacs/options/ioptionscontainer.h"
  60 #include "gromacs/options/timeunitmanager.h"
  61 #include "gromacs/selection/selectioncollection.h"
  62 #include "gromacs/utility/exceptions.h"
  63 #include "gromacs/utility/gmxassert.h"
  64 #include "gromacs/utility/programcontext.h"
  65 #include "gromacs/utility/stringutil.h"
  66 #include "gromacs/utility/unique_cptr.h"
  67
  68 namespace gmx
  69 {
  70
  71 /********************************************************************
  72  * AnalysisDataPlotSettings
  73  */
  74
  75 AnalysisDataPlotSettings::AnalysisDataPlotSettings() :
  76     selections_(nullptr), timeUnit_(TimeUnit::Default), plotFormat_(XvgFormat::Xmgrace)
  77 {
  78 }
  79
  80 void AnalysisDataPlotSettings::setSelectionCollection(const SelectionCollection* selections)
  81 {
  82     selections_ = selections;
  83 }
  84
  85 /*! \brief Names for XvgFormat
  86  *
  87  * Technically this duplicates a definition in pargs.cpp for legacy
  88  * support code, but as the latter will go away and the alternatives
  89  * are ugly, the duplication is acceptable. */
  90 const gmx::EnumerationArray<XvgFormat, const char*> c_xvgFormatNames = {
  91     { "xmgrace", "xmgr", "none" }
  92 };
  93
  94 void AnalysisDataPlotSettings::initOptions(IOptionsContainer* options)
  95 {
  96     options->addOption(
  97             EnumOption<XvgFormat>("xvg").enumValue(c_xvgFormatNames).store(&plotFormat_).description("Plot formatting"));
  98 }
  99
 100
 101 /********************************************************************
 102  * AbstractPlotModule::Impl
 103  */
 104
 105 class AbstractPlotModule::Impl
 106 {
 107 public:
 108     explicit Impl(const AnalysisDataPlotSettings& settings);
 109     ~Impl();
 110
 111     void closeFile();
 112
 113     AnalysisDataPlotSettings settings_;
 114     std::string              filename_;
 115     FILE*                    fp_;
 116
 117     bool                     bPlain_;
 118     bool                     bOmitX_;
 119     bool                     bErrorsAsSeparateColumn_;
 120     std::string              title_;
 121     std::string              subtitle_;
 122     std::string              xlabel_;
 123     std::string              ylabel_;
 124     std::vector<std::string> legend_;
 125     std::string              xformat_;
 126     std::string              yformat_;
 127     real                     xscale_;
 128 };
 129
 130 AbstractPlotModule::Impl::Impl(const AnalysisDataPlotSettings& settings) :
 131     settings_(settings),
 132     fp_(nullptr),
 133     bPlain_(false),
 134     bOmitX_(false),
 135     bErrorsAsSeparateColumn_(false),
 136     xformat_("%11.3f"),
 137     yformat_(" %8.3f"),
 138     xscale_(1.0)
 139 {
 140 }
 141
 142 AbstractPlotModule::Impl::~Impl()
 143 {
 144     closeFile();
 145 }
 146
 147
 148 void AbstractPlotModule::Impl::closeFile()
 149 {
 150     if (fp_ != nullptr)
 151     {
 152         if (bPlain_)
 153         {
 154             gmx_fio_fclose(fp_);
 155         }
 156         else
 157         {
 158             xvgrclose(fp_);
 159         }
 160         fp_ = nullptr;
 161     }
 162 }
 163
 164
 165 /********************************************************************
 166  * AbstractPlotModule
 167  */
 168 /*! \cond libapi */
 169 AbstractPlotModule::AbstractPlotModule() : impl_(new Impl(AnalysisDataPlotSettings())) {}
 170
 171 AbstractPlotModule::AbstractPlotModule(const AnalysisDataPlotSettings& settings) :
 172     impl_(new Impl(settings))
 173 {
 174 }
 175 //! \endcond
 176
 177 AbstractPlotModule::~AbstractPlotModule() {}
 178
 179
 180 void AbstractPlotModule::setSettings(const AnalysisDataPlotSettings& settings)
 181 {
 182     impl_->settings_ = settings;
 183 }
 184
 185
 186 void AbstractPlotModule::setFileName(const std::string& filename)
 187 {
 188     impl_->filename_ = filename;
 189 }
 190
 191
 192 void AbstractPlotModule::setPlainOutput(bool bPlain)
 193 {
 194     impl_->bPlain_ = bPlain;
 195 }
 196
 197
 198 void AbstractPlotModule::setErrorsAsSeparateColumn(bool bSeparate)
 199 {
 200     impl_->bErrorsAsSeparateColumn_ = bSeparate;
 201 }
 202
 203
 204 void AbstractPlotModule::setOmitX(bool bOmitX)
 205 {
 206     impl_->bOmitX_ = bOmitX;
 207 }
 208
 209
 210 void AbstractPlotModule::setTitle(const char* title)
 211 {
 212     impl_->title_ = title;
 213 }
 214
 215 void AbstractPlotModule::setTitle(const std::string& title)
 216 {
 217     impl_->title_ = title;
 218 }
 219
 220
 221 void AbstractPlotModule::setSubtitle(const char* subtitle)
 222 {
 223     impl_->subtitle_ = subtitle;
 224 }
 225
 226
 227 void AbstractPlotModule::setSubtitle(const std::string& subtitle)
 228 {
 229     impl_->subtitle_ = subtitle;
 230 }
 231
 232
 233 void AbstractPlotModule::setXLabel(const char* label)
 234 {
 235     impl_->xlabel_ = label;
 236 }
 237
 238
 239 void AbstractPlotModule::setXAxisIsTime()
 240 {
 241     TimeUnitManager manager(impl_->settings_.timeUnit());
 242     impl_->xlabel_ = formatString("Time (%s)", manager.timeUnitAsString());
 243     impl_->xscale_ = manager.inverseTimeScaleFactor();
 244 }
 245
 246
 247 void AbstractPlotModule::setYLabel(const char* label)
 248 {
 249     impl_->ylabel_ = label;
 250 }
 251
 252
 253 void AbstractPlotModule::setLegend(int nsets, const char* const* setname)
 254 {
 255     impl_->legend_.reserve(impl_->legend_.size() + nsets);
 256     for (int i = 0; i < nsets; ++i)
 257     {
 258         appendLegend(setname[i]);
 259     }
 260 }
 261
 262
 263 void AbstractPlotModule::appendLegend(const char* setname)
 264 {
 265     impl_->legend_.emplace_back(setname);
 266 }
 267
 268
 269 void AbstractPlotModule::appendLegend(const std::string& setname)
 270 {
 271     impl_->legend_.push_back(setname);
 272 }
 273
 274
 275 void AbstractPlotModule::setXFormat(int width, int precision, char format)
 276 {
 277     GMX_RELEASE_ASSERT(width >= 0 && precision >= 0 && width <= 99 && precision <= 99,
 278                        "Invalid width or precision");
 279     GMX_RELEASE_ASSERT(strchr("eEfFgG", format) != nullptr, "Invalid format specifier");
 280     impl_->xformat_ = formatString("%%%d.%d%c", width, precision, format);
 281 }
 282
 283
 284 void AbstractPlotModule::setYFormat(int width, int precision, char format)
 285 {
 286     GMX_RELEASE_ASSERT(width >= 0 && precision >= 0 && width <= 99 && precision <= 99,
 287                        "Invalid width or precision");
 288     GMX_RELEASE_ASSERT(strchr("eEfFgG", format) != nullptr, "Invalid format specifier");
 289     impl_->yformat_ = formatString(" %%%d.%d%c", width, precision, format);
 290 }
 291
 292
 293 int AbstractPlotModule::flags() const
 294 {
 295     return efAllowMissing | efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 296 }
 297
 298
 299 void AbstractPlotModule::dataStarted(AbstractAnalysisData* /* data */)
 300 {
 301     if (!impl_->filename_.empty())
 302     {
 303         if (impl_->bPlain_)
 304         {
 305             impl_->fp_ = gmx_fio_fopen(impl_->filename_.c_str(), "w");
 306         }
 307         else
 308         {
 309             const TimeUnit    timeUnit  = impl_->settings_.timeUnit();
 310             const XvgFormat   xvgFormat = impl_->settings_.plotFormat();
 311             gmx_output_env_t* oenv;
 312             output_env_init(&oenv, getProgramContext(), timeUnit, FALSE, xvgFormat, 0);
 313             const unique_cptr<gmx_output_env_t, output_env_done> oenvGuard(oenv);
 314             impl_->fp_ = xvgropen(
 315                     impl_->filename_.c_str(), impl_->title_.c_str(), impl_->xlabel_, impl_->ylabel_, oenv);
 316             const SelectionCollection* selections = impl_->settings_.selectionCollection();
 317             if (selections != nullptr && output_env_get_xvg_format(oenv) != XvgFormat::None)
 318             {
 319                 selections->printXvgrInfo(impl_->fp_);
 320             }
 321             if (!impl_->subtitle_.empty())
 322             {
 323                 xvgr_subtitle(impl_->fp_, impl_->subtitle_.c_str(), oenv);
 324             }
 325             if (output_env_get_print_xvgr_codes(oenv) && !impl_->legend_.empty())
 326             {
 327                 std::vector<const char*> legend;
 328                 legend.reserve(impl_->legend_.size());
 329                 for (size_t i = 0; i < impl_->legend_.size(); ++i)
 330                 {
 331                     legend.push_back(impl_->legend_[i].c_str());
 332                 }
 333                 xvgr_legend(impl_->fp_, legend.size(), legend.data(), oenv);
 334             }
 335         }
 336     }
 337 }
 338
 339
 340 void AbstractPlotModule::frameStarted(const AnalysisDataFrameHeader& header)
 341 {
 342     if (!isFileOpen())
 343     {
 344         return;
 345     }
 346     if (!impl_->bOmitX_)
 347     {
 348         std::fprintf(impl_->fp_, impl_->xformat_.c_str(), header.x() * impl_->xscale_);
 349     }
 350 }
 351
 352
 353 void AbstractPlotModule::frameFinished(const AnalysisDataFrameHeader& /*header*/)
 354 {
 355     if (!isFileOpen())
 356     {
 357         return;
 358     }
 359     std::fprintf(impl_->fp_, "\n");
 360 }
 361
 362
 363 void AbstractPlotModule::dataFinished()
 364 {
 365     impl_->closeFile();
 366 }
 367
 368 /*! \cond libapi */
 369 bool AbstractPlotModule::isFileOpen() const
 370 {
 371     return impl_->fp_ != nullptr;
 372 }
 373
 374
 375 void AbstractPlotModule::writeValue(const AnalysisDataValue& value) const
 376 {
 377     GMX_ASSERT(isFileOpen(), "File not opened, but write attempted");
 378     const real y = value.isSet() ? value.value() : 0.0;
 379     std::fprintf(impl_->fp_, impl_->yformat_.c_str(), y);
 380     if (impl_->bErrorsAsSeparateColumn_)
 381     {
 382         const real dy = value.isSet() ? value.error() : 0.0;
 383         std::fprintf(impl_->fp_, impl_->yformat_.c_str(), dy);
 384     }
 385 }
 386 //! \endcond
 387
 388 /********************************************************************
 389  * DataPlotModule
 390  */
 391
 392 AnalysisDataPlotModule::AnalysisDataPlotModule() {}
 393
 394 AnalysisDataPlotModule::AnalysisDataPlotModule(const AnalysisDataPlotSettings& settings) :
 395     AbstractPlotModule(settings)
 396 {
 397 }
 398
 399
 400 void AnalysisDataPlotModule::pointsAdded(const AnalysisDataPointSetRef& points)
 401 {
 402     if (!isFileOpen())
 403     {
 404         return;
 405     }
 406     for (int i = 0; i < points.columnCount(); ++i)
 407     {
 408         writeValue(points.values()[i]);
 409     }
 410 }
 411
 412
 413 /********************************************************************
 414  * DataVectorPlotModule
 415  */
 416
 417 AnalysisDataVectorPlotModule::AnalysisDataVectorPlotModule() : bWrite_{ true, true, true, false } {}
 418
 419
 420 AnalysisDataVectorPlotModule::AnalysisDataVectorPlotModule(const AnalysisDataPlotSettings& settings) :
 421     AbstractPlotModule(settings), bWrite_{ true, true, true, false }
 422 {
 423 }
 424
 425
 426 void AnalysisDataVectorPlotModule::setWriteX(bool bWrite)
 427 {
 428     bWrite_[XX] = bWrite;
 429 }
 430
 431
 432 void AnalysisDataVectorPlotModule::setWriteY(bool bWrite)
 433 {
 434     bWrite_[YY] = bWrite;
 435 }
 436
 437
 438 void AnalysisDataVectorPlotModule::setWriteZ(bool bWrite)
 439 {
 440     bWrite_[ZZ] = bWrite;
 441 }
 442
 443
 444 void AnalysisDataVectorPlotModule::setWriteNorm(bool bWrite)
 445 {
 446     bWrite_[DIM] = bWrite;
 447 }
 448
 449
 450 void AnalysisDataVectorPlotModule::setWriteMask(const bool bWrite[DIM + 1])
 451 {
 452     for (int i = 0; i < DIM + 1; ++i)
 453     {
 454         bWrite_[i] = bWrite[i];
 455     }
 456 }
 457
 458
 459 void AnalysisDataVectorPlotModule::pointsAdded(const AnalysisDataPointSetRef& points)
 460 {
 461     if (points.firstColumn() % DIM != 0 || points.columnCount() % DIM != 0)
 462     {
 463         GMX_THROW(APIError("Partial data points"));
 464     }
 465     if (!isFileOpen())
 466     {
 467         return;
 468     }
 469     for (int i = 0; i < points.columnCount(); i += 3)
 470     {
 471         for (int d = 0; d < DIM; ++d)
 472         {
 473             if (bWrite_[d])
 474             {
 475                 writeValue(points.values()[i + d]);
 476             }
 477         }
 478         if (bWrite_[DIM])
 479         {
 480             const rvec        y = { points.y(i), points.y(i + 1), points.y(i + 2) };
 481             AnalysisDataValue value(norm(y));
 482             writeValue(value);
 483         }
 484     }
 485 }
 486
 487 } // namespace gmx