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