src/gromacs/analysisdata/modules/plot.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2011,2012,2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements classes in plot.h.
  38  *
  39  * \ingroup module_analysisdata
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  */
  42 #include "gromacs/analysisdata/modules/plot.h"
  43
  44 #include <string>
  45 #include <vector>
  46
  47 #include <cstdio>
  48 #include <cstring>
  49
  50 #include <boost/shared_ptr.hpp>
  51
  52 #include "gromacs/legacyheaders/gmxfio.h"
  53 #include "gromacs/legacyheaders/oenv.h"
  54 #include "gromacs/legacyheaders/vec.h"
  55 #include "gromacs/legacyheaders/xvgr.h"
  56
  57 #include "gromacs/options/basicoptions.h"
  58 #include "gromacs/analysisdata/dataframe.h"
  59 #include "gromacs/options/options.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/stringutil.h"
  65
  66 namespace
  67 {
  68
  69 //! Enum values for plot formats.
  70 const char *const g_plotFormats[] = {
  71     "none", "xmgrace", "xmgr"
  72 };
  73
  74 } // namespace
  75
  76 namespace gmx
  77 {
  78
  79 /********************************************************************
  80  * AnalysisDataPlotSettings
  81  */
  82
  83 AnalysisDataPlotSettings::AnalysisDataPlotSettings()
  84     : selections_(NULL), timeUnit_(eTimeUnit_ps), plotFormat_(1)
  85 {
  86 }
  87
  88 void
  89 AnalysisDataPlotSettings::setSelectionCollection(const SelectionCollection *selections)
  90 {
  91     selections_ = selections;
  92 }
  93
  94
  95 void
  96 AnalysisDataPlotSettings::addOptions(Options *options)
  97 {
  98     options->addOption(StringOption("xvg").enumValue(g_plotFormats)
  99                            .defaultValue("xmgrace")
 100                            .storeEnumIndex(&plotFormat_)
 101                            .description("Plot formatting"));
 102 }
 103
 104
 105 /********************************************************************
 106  * AbstractPlotModule::Impl
 107  */
 108
 109 class AbstractPlotModule::Impl
 110 {
 111     public:
 112         explicit Impl(const AnalysisDataPlotSettings &settings);
 113         ~Impl();
 114
 115         void closeFile();
 116
 117         AnalysisDataPlotSettings  settings_;
 118         std::string               filename_;
 119         FILE                     *fp_;
 120
 121         bool                      bPlain_;
 122         bool                      gOmitX_;
 123         std::string               title_;
 124         std::string               subtitle_;
 125         std::string               xlabel_;
 126         std::string               ylabel_;
 127         std::vector<std::string>  legend_;
 128         char                      xformat_[15];
 129         char                      yformat_[15];
 130         real                      xscale_;
 131 };
 132
 133 AbstractPlotModule::Impl::Impl(const AnalysisDataPlotSettings &settings)
 134     : settings_(settings), fp_(NULL), bPlain_(false), gOmitX_(false),
 135       xscale_(1.0)
 136 {
 137     strcpy(xformat_, "%11.3f");
 138     strcpy(yformat_, " %8.3f");
 139 }
 140
 141 AbstractPlotModule::Impl::~Impl()
 142 {
 143     closeFile();
 144 }
 145
 146
 147 void
 148 AbstractPlotModule::Impl::closeFile()
 149 {
 150     if (fp_ != NULL)
 151     {
 152         if (bPlain_)
 153         {
 154             gmx_fio_fclose(fp_);
 155         }
 156         else
 157         {
 158             xvgrclose(fp_);
 159         }
 160         fp_ = NULL;
 161     }
 162 }
 163
 164
 165 /********************************************************************
 166  * AbstractPlotModule
 167  */
 168 /*! \cond libapi */
 169 AbstractPlotModule::AbstractPlotModule()
 170     : impl_(new Impl(AnalysisDataPlotSettings()))
 171 {
 172 }
 173
 174 AbstractPlotModule::AbstractPlotModule(const AnalysisDataPlotSettings &settings)
 175     : impl_(new Impl(settings))
 176 {
 177 }
 178 //! \endcond
 179
 180 AbstractPlotModule::~AbstractPlotModule()
 181 {
 182 }
 183
 184
 185 void
 186 AbstractPlotModule::setSettings(const AnalysisDataPlotSettings &settings)
 187 {
 188     impl_->settings_ = settings;
 189 }
 190
 191
 192 void
 193 AbstractPlotModule::setFileName(const std::string &filename)
 194 {
 195     impl_->filename_ = filename;
 196 }
 197
 198
 199 void
 200 AbstractPlotModule::setPlainOutput(bool bPlain)
 201 {
 202     impl_->bPlain_ = bPlain;
 203 }
 204
 205
 206 void
 207 AbstractPlotModule::setOmitX(bool bOmitX)
 208 {
 209     impl_->gOmitX_ = bOmitX;
 210 }
 211
 212
 213 void
 214 AbstractPlotModule::setTitle(const char *title)
 215 {
 216     impl_->title_ = title;
 217 }
 218
 219
 220 void
 221 AbstractPlotModule::setSubtitle(const char *subtitle)
 222 {
 223     impl_->subtitle_ = subtitle;
 224 }
 225
 226
 227 void
 228 AbstractPlotModule::setXLabel(const char *label)
 229 {
 230     impl_->xlabel_ = label;
 231 }
 232
 233
 234 void
 235 AbstractPlotModule::setXAxisIsTime()
 236 {
 237     TimeUnitManager manager(impl_->settings_.timeUnit());
 238     impl_->xlabel_ = formatString("Time (%s)", manager.timeUnitAsString());
 239     impl_->xscale_ = manager.inverseTimeScaleFactor();
 240 }
 241
 242
 243 void
 244 AbstractPlotModule::setYLabel(const char *label)
 245 {
 246     impl_->ylabel_ = label;
 247 }
 248
 249
 250 void
 251 AbstractPlotModule::setLegend(int nsets, const char * const *setname)
 252 {
 253     impl_->legend_.reserve(impl_->legend_.size() + nsets);
 254     for (int i = 0; i < nsets; ++i)
 255     {
 256         appendLegend(setname[i]);
 257     }
 258 }
 259
 260
 261 void
 262 AbstractPlotModule::appendLegend(const char *setname)
 263 {
 264     impl_->legend_.push_back(setname);
 265 }
 266
 267
 268 void
 269 AbstractPlotModule::setXFormat(int width, int precision, char format)
 270 {
 271     GMX_RELEASE_ASSERT(width >= 0 && precision >= 0
 272                        && width <= 99 && precision <= 99,
 273                        "Invalid width or precision");
 274     GMX_RELEASE_ASSERT(strchr("eEfFgG", format) != NULL,
 275                        "Invalid format specifier");
 276     sprintf(impl_->xformat_, "%%%d.%d%c", width, precision, format);
 277 }
 278
 279
 280 void
 281 AbstractPlotModule::setYFormat(int width, int precision, char format)
 282 {
 283     GMX_RELEASE_ASSERT(width >= 0 && precision >= 0
 284                        && width <= 99 && precision <= 99,
 285                        "Invalid width or precision");
 286     GMX_RELEASE_ASSERT(strchr("eEfFgG", format) != NULL,
 287                        "Invalid format specifier");
 288     sprintf(impl_->yformat_, " %%%d.%d%c", width, precision, format);
 289 }
 290
 291
 292 int
 293 AbstractPlotModule::flags() const
 294 {
 295     return efAllowMulticolumn | efAllowMultipoint;
 296 }
 297
 298
 299 void
 300 AbstractPlotModule::dataStarted(AbstractAnalysisData *data)
 301 {
 302     if (!impl_->filename_.empty())
 303     {
 304         if (impl_->bPlain_)
 305         {
 306             impl_->fp_ = gmx_fio_fopen(impl_->filename_.c_str(), "w");
 307         }
 308         else
 309         {
 310             time_unit_t  time_unit
 311                 = static_cast<time_unit_t>(impl_->settings_.timeUnit() + 1);
 312             xvg_format_t xvg_format
 313                 = (impl_->settings_.plotFormat() > 0
 314                    ? static_cast<xvg_format_t>(impl_->settings_.plotFormat())
 315                    : exvgNONE);
 316             output_env_t                  oenv;
 317             output_env_init(&oenv, 0, NULL, time_unit, FALSE, xvg_format, 0, 0);
 318             boost::shared_ptr<output_env> oenvGuard(oenv, &output_env_done);
 319             impl_->fp_ = xvgropen(impl_->filename_.c_str(), impl_->title_.c_str(),
 320                                   impl_->xlabel_.c_str(), impl_->ylabel_.c_str(),
 321                                   oenv);
 322             const SelectionCollection *selections
 323                 = impl_->settings_.selectionCollection();
 324             if (selections != NULL)
 325             {
 326                 selections->printXvgrInfo(impl_->fp_, oenv);
 327             }
 328             if (!impl_->subtitle_.empty())
 329             {
 330                 xvgr_subtitle(impl_->fp_, impl_->subtitle_.c_str(), oenv);
 331             }
 332             if (output_env_get_print_xvgr_codes(oenv)
 333                 && !impl_->legend_.empty())
 334             {
 335                 std::vector<const char *> legend;
 336                 legend.reserve(impl_->legend_.size());
 337                 for (size_t i = 0; i < impl_->legend_.size(); ++i)
 338                 {
 339                     legend.push_back(impl_->legend_[i].c_str());
 340                 }
 341                 xvgr_legend(impl_->fp_, legend.size(), &legend[0], oenv);
 342             }
 343         }
 344     }
 345 }
 346
 347
 348 void
 349 AbstractPlotModule::frameStarted(const AnalysisDataFrameHeader &frame)
 350 {
 351     if (!isFileOpen())
 352     {
 353         return;
 354     }
 355     if (!impl_->gOmitX_)
 356     {
 357         std::fprintf(impl_->fp_, impl_->xformat_, frame.x() * impl_->xscale_);
 358     }
 359 }
 360
 361
 362 void
 363 AbstractPlotModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 364 {
 365     if (!isFileOpen())
 366     {
 367         return;
 368     }
 369     std::fprintf(impl_->fp_, "\n");
 370 }
 371
 372
 373 void
 374 AbstractPlotModule::dataFinished()
 375 {
 376     impl_->closeFile();
 377 }
 378
 379 /*! \cond libapi */
 380 bool
 381 AbstractPlotModule::isFileOpen() const
 382 {
 383     return impl_->fp_ != NULL;
 384 }
 385
 386
 387 void
 388 AbstractPlotModule::writeValue(real value) const
 389 {
 390     GMX_ASSERT(isFileOpen(), "File not opened, but write attempted");
 391     std::fprintf(impl_->fp_, impl_->yformat_, value);
 392 }
 393 //! \endcond
 394
 395 /********************************************************************
 396  * DataPlotModule
 397  */
 398
 399 AnalysisDataPlotModule::AnalysisDataPlotModule()
 400 {
 401 }
 402
 403 AnalysisDataPlotModule::AnalysisDataPlotModule(
 404         const AnalysisDataPlotSettings &settings)
 405     : AbstractPlotModule(settings)
 406 {
 407 }
 408
 409
 410 void
 411 AnalysisDataPlotModule::pointsAdded(const AnalysisDataPointSetRef &points)
 412 {
 413     if (!isFileOpen())
 414     {
 415         return;
 416     }
 417     for (int i = 0; i < points.columnCount(); ++i)
 418     {
 419         writeValue(points.y(i));
 420     }
 421 }
 422
 423
 424 /********************************************************************
 425  * DataVectorPlotModule
 426  */
 427
 428 AnalysisDataVectorPlotModule::AnalysisDataVectorPlotModule()
 429 {
 430     for (int i = 0; i < DIM; ++i)
 431     {
 432         bWrite_[i] = true;
 433     }
 434     bWrite_[DIM] = false;
 435 }
 436
 437
 438 AnalysisDataVectorPlotModule::AnalysisDataVectorPlotModule(
 439         const AnalysisDataPlotSettings &settings)
 440     : AbstractPlotModule(settings)
 441 {
 442     for (int i = 0; i < DIM; ++i)
 443     {
 444         bWrite_[i] = true;
 445     }
 446     bWrite_[DIM] = false;
 447 }
 448
 449
 450 void
 451 AnalysisDataVectorPlotModule::setWriteX(bool bWrite)
 452 {
 453     bWrite_[XX] = bWrite;
 454 }
 455
 456
 457 void
 458 AnalysisDataVectorPlotModule::setWriteY(bool bWrite)
 459 {
 460     bWrite_[YY] = bWrite;
 461 }
 462
 463
 464 void
 465 AnalysisDataVectorPlotModule::setWriteZ(bool bWrite)
 466 {
 467     bWrite_[ZZ] = bWrite;
 468 }
 469
 470
 471 void
 472 AnalysisDataVectorPlotModule::setWriteNorm(bool bWrite)
 473 {
 474     bWrite_[DIM] = bWrite;
 475 }
 476
 477
 478 void
 479 AnalysisDataVectorPlotModule::setWriteMask(bool bWrite[DIM + 1])
 480 {
 481     for (int i = 0; i < DIM + 1; ++i)
 482     {
 483         bWrite_[i] = bWrite[i];
 484     }
 485 }
 486
 487
 488 void
 489 AnalysisDataVectorPlotModule::pointsAdded(const AnalysisDataPointSetRef &points)
 490 {
 491     if (points.firstColumn() % DIM != 0)
 492     {
 493         GMX_THROW(APIError("Partial data points"));
 494     }
 495     if (!isFileOpen())
 496     {
 497         return;
 498     }
 499     for (int i = 0; i < points.columnCount(); i += 3)
 500     {
 501         for (int d = 0; d < DIM; ++d)
 502         {
 503             if (bWrite_[i])
 504             {
 505                 writeValue(points.y(i + d));
 506             }
 507         }
 508         if (bWrite_[DIM])
 509         {
 510             rvec y = { points.y(i), points.y(i + 1), points.y(i + 2) };
 511             writeValue(norm(y));
 512         }
 513     }
 514 }
 515
 516 } // namespace gmx