This provides a lot more flexible alternative to g_dist -lt.
Added an analysisdata module that computes lifetime statistics and test
for it. Currently tuned for this one application, but should be
possible to make more flexible if required. Using this module, the
implementation in gmx-select is straightforward.
Part of #665.
Change-Id: I8bff89599c9d8a820c761b426c01ceafa9afcc5a
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2013, by the GROMACS development team, led by
+ * David van der Spoel, Berk Hess, Erik Lindahl, and including many
+ * others, as listed in the AUTHORS file in the top-level source
+ * directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements gmx::AnalysisDataLifetimeModule.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_analysisdata
+ */
+#include "lifetime.h"
+
+#include <cmath>
+
+#include <deque>
+#include <vector>
+
+#include "gromacs/analysisdata/dataframe.h"
+#include "gromacs/analysisdata/datastorage.h"
+
+namespace gmx
+{
+
+/********************************************************************
+ * AnalysisDataLifetimeModule
+ */
+
+/*! \internal \brief
+ * Private implementation class for AnalysisDataLifetimeModule.
+ *
+ * \ingroup module_analysisdata
+ */
+class AnalysisDataLifetimeModule::Impl
+{
+ public:
+ //! Container type for storing a histogram during the calculation.
+ typedef std::deque<int> LifetimeHistogram;
+
+ //! Initializes the implementation class with empty/default values.
+ Impl() : firstx_(0.0), lastx_(0.0), frameCount_(0), bCumulative_(false)
+ {
+ }
+
+ /*! \brief
+ * Increments a lifetime histogram with a single lifetime.
+ *
+ * \param[in] dataSet Index of the histogram to increment.
+ * \param[in] lifetime Lifetime to add to the histogram.
+ */
+ void addLifetime(int dataSet, int lifetime)
+ {
+ if (lifetime > 0)
+ {
+ LifetimeHistogram &histogram = lifetimeHistograms_[dataSet];
+ if (histogram.size() < static_cast<unsigned>(lifetime))
+ {
+ histogram.resize(lifetime, 0);
+ }
+ ++histogram[lifetime - 1];
+ }
+ }
+
+ //! X value of the first frame (used for determining output spacing).
+ real firstx_;
+ //! X value of the last frame (used for determining output spacing).
+ real lastx_;
+ //! Total number of frames (used for normalization and output spacing).
+ int frameCount_;
+ //! Whether to add subintervals of longer intervals explicitly.
+ bool bCumulative_;
+ /*! \brief
+ * Length of current continuously present interval for each data column.
+ *
+ * While frame N has been processed, stores the length of an interval
+ * for each data column where that column has been continuously present
+ * up to and including frame N.
+ */
+ std::vector<std::vector<int> > currentLifetimes_;
+ /*! \brief
+ * Accumulated lifetime histograms for each data set.
+ */
+ std::vector<LifetimeHistogram> lifetimeHistograms_;
+};
+
+AnalysisDataLifetimeModule::AnalysisDataLifetimeModule()
+ : impl_(new Impl())
+{
+}
+
+AnalysisDataLifetimeModule::~AnalysisDataLifetimeModule()
+{
+}
+
+void AnalysisDataLifetimeModule::setCumulative(bool bCumulative)
+{
+ impl_->bCumulative_ = bCumulative;
+}
+
+int AnalysisDataLifetimeModule::flags() const
+{
+ return efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
+}
+
+void
+AnalysisDataLifetimeModule::dataStarted(AbstractAnalysisData *data)
+{
+ impl_->currentLifetimes_.reserve(data->dataSetCount());
+ impl_->lifetimeHistograms_.reserve(data->dataSetCount());
+ for (int i = 0; i < data->dataSetCount(); ++i)
+ {
+ impl_->currentLifetimes_.push_back(std::vector<int>(data->columnCount(i), 0));
+ impl_->lifetimeHistograms_.push_back(std::deque<int>());
+ }
+}
+
+void
+AnalysisDataLifetimeModule::frameStarted(const AnalysisDataFrameHeader &header)
+{
+ if (header.index() == 0)
+ {
+ impl_->firstx_ = header.x();
+ }
+ impl_->lastx_ = header.x();
+ ++impl_->frameCount_;
+ // TODO: Check the input for even spacing.
+}
+
+void
+AnalysisDataLifetimeModule::pointsAdded(const AnalysisDataPointSetRef &points)
+{
+ const int dataSet = points.dataSetIndex();
+ // This assumption is strictly not necessary, but this is how the
+ // framework works currently, and makes the code below simpler.
+ GMX_ASSERT(points.firstColumn() == 0
+ && points.lastColumn() == static_cast<int>(impl_->currentLifetimes_[dataSet].size()) - 1,
+ "Point set should cover all columns");
+ for (int i = 0; i < points.columnCount(); ++i)
+ {
+ // TODO: Perhaps add control over how this is determined?
+ const bool bPresent = points.present(i) && points.y(i) > 0.0;
+ if (bPresent)
+ {
+ ++impl_->currentLifetimes_[dataSet][i];
+ }
+ else if (impl_->currentLifetimes_[dataSet][i] > 0)
+ {
+ impl_->addLifetime(dataSet, impl_->currentLifetimes_[dataSet][i]);
+ impl_->currentLifetimes_[dataSet][i] = 0;
+ }
+ }
+}
+
+void
+AnalysisDataLifetimeModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
+{
+}
+
+void
+AnalysisDataLifetimeModule::dataFinished()
+{
+ // Need to process the elements present in the last frame explicitly.
+ for (size_t i = 0; i < impl_->currentLifetimes_.size(); ++i)
+ {
+ for (size_t j = 0; j < impl_->currentLifetimes_[i].size(); ++j)
+ {
+ impl_->addLifetime(i, impl_->currentLifetimes_[i][j]);
+ }
+ }
+ impl_->currentLifetimes_.clear();
+
+ if (impl_->bCumulative_)
+ {
+ // Sum up subintervals of longer intervals into the histograms
+ // if explicitly requested.
+ std::vector<Impl::LifetimeHistogram>::iterator histogram;
+ for (histogram = impl_->lifetimeHistograms_.begin();
+ histogram != impl_->lifetimeHistograms_.end();
+ ++histogram)
+ {
+ Impl::LifetimeHistogram::iterator shorter, longer;
+ for (shorter = histogram->begin(); shorter != histogram->end(); ++shorter)
+ {
+ int subIntervalCount = 2;
+ for (longer = shorter + 1; longer != histogram->end();
+ ++longer, ++subIntervalCount)
+ {
+ // Interval of length shorter contains (longer - shorter + 1)
+ // continuous intervals of length longer.
+ *shorter += subIntervalCount * (*longer);
+ }
+ }
+ }
+ }
+
+ // X spacing is determined by averaging from the first and last frame
+ // instead of first two frames to avoid rounding issues.
+ const real spacing =
+ (impl_->frameCount_ > 1)
+ ? (impl_->lastx_ - impl_->firstx_) / (impl_->frameCount_ - 1)
+ : 0.0;
+ setXAxis(0.0, spacing);
+
+ // Determine output dimensionality to cover all the histograms.
+ setColumnCount(impl_->lifetimeHistograms_.size());
+ std::vector<Impl::LifetimeHistogram>::const_iterator histogram;
+ size_t maxLifetime = 1;
+ for (histogram = impl_->lifetimeHistograms_.begin();
+ histogram != impl_->lifetimeHistograms_.end();
+ ++histogram)
+ {
+ maxLifetime = std::max(maxLifetime, histogram->size());
+ }
+ setRowCount(maxLifetime);
+
+ // Fill up the output data from the histograms.
+ allocateValues();
+ int column = 0;
+ for (histogram = impl_->lifetimeHistograms_.begin();
+ histogram != impl_->lifetimeHistograms_.end();
+ ++histogram, ++column)
+ {
+ int row = 0;
+ Impl::LifetimeHistogram::const_iterator i;
+ for (i = histogram->begin(); i != histogram->end(); ++i, ++row)
+ {
+ // Normalize by the number of frames, taking into account the
+ // length of the interval (interval of length N cannot start in
+ // N-1 last frames). row is always smaller than frameCount_
+ // because of the histograms have at most frameCount_ entries.
+ const real normalized = *i / static_cast<real>(impl_->frameCount_ - row);
+ value(row, column).setValue(normalized);
+ }
+ // Pad the rest of the histogram with zeros to match the longest
+ // histogram.
+ for (; row < rowCount(); ++row)
+ {
+ value(row, column).setValue(0.0);
+ }
+ }
+ impl_->lifetimeHistograms_.clear();
+ valuesReady();
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2013, by the GROMACS development team, led by
+ * David van der Spoel, Berk Hess, Erik Lindahl, and including many
+ * others, as listed in the AUTHORS file in the top-level source
+ * directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \file
+ * \brief
+ * Declares gmx::AnalysisDataLifetimeModule.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \inpublicapi
+ * \ingroup module_analysisdata
+ */
+#ifndef GMX_ANALYSISDATA_MODULES_LIFETIME_H
+#define GMX_ANALYSISDATA_MODULES_LIFETIME_H
+
+#include "../arraydata.h"
+#include "../datamodule.h"
+#include "../../utility/common.h"
+
+namespace gmx
+{
+
+/*! \brief
+ * Data module for computing lifetime histograms for columns in input data.
+ *
+ * The input data set is treated as a boolean array: each value that is present
+ * (AnalysisDataValue::isPresent() returns true) and is >0 is treated as
+ * present, other values are treated as absent.
+ * For each input data set, analyzes the columns to identify the intervals
+ * where a column is continuously present.
+ * Produces a histogram from the lengths of these intervals.
+ * Input data should have frames with evenly spaced x values.
+ *
+ * Output data contains one column for each data set in the input data.
+ * This column gives the lifetime histogram for the corresponding data set.
+ * x axis in the output is spaced the same as in the input data, and extends
+ * as long as required to cover all the histograms.
+ * Histograms are padded with zeros as required to be of the same length.
+ * setCumulative() can be used to alter the handling of subintervals in the
+ * output histogram.
+ *
+ * The output data becomes available only after the input data has been
+ * finished.
+ *
+ * \inpublicapi
+ * \ingroup module_analysisdata
+ */
+class AnalysisDataLifetimeModule : public AbstractAnalysisArrayData,
+ public AnalysisDataModuleInterface
+{
+ public:
+ AnalysisDataLifetimeModule();
+ virtual ~AnalysisDataLifetimeModule();
+
+ /*! \brief
+ * Sets a cumulative histogram mode.
+ *
+ * \param[in] bCumulative If true, all subintervals of a long
+ * interval are also explicitly added into the histogram.
+ *
+ * Does not throw.
+ */
+ void setCumulative(bool bCumulative);
+
+ virtual int flags() const;
+
+ virtual void dataStarted(AbstractAnalysisData *data);
+ virtual void frameStarted(const AnalysisDataFrameHeader &header);
+ virtual void pointsAdded(const AnalysisDataPointSetRef &points);
+ virtual void frameFinished(const AnalysisDataFrameHeader &header);
+ virtual void dataFinished();
+
+ private:
+ class Impl;
+
+ PrivateImplPointer<Impl> impl_;
+};
+
+//! Smart pointer to manage an AnalysisDataLifetimeModule object.
+typedef boost::shared_ptr<AnalysisDataLifetimeModule>
+ AnalysisDataLifetimeModulePointer;
+
+} // namespace gmx
+
+#endif
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2011,2012, by the GROMACS development team, led by
+# Copyright (c) 2011,2012,2013, by the GROMACS development team, led by
# David van der Spoel, Berk Hess, Erik Lindahl, and including many
# others, as listed in the AUTHORS file in the top-level source
# directory and at http://www.gromacs.org.
analysisdata.cpp
arraydata.cpp
average.cpp
- histogram.cpp)
+ histogram.cpp
+ lifetime.cpp)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2013, by the GROMACS development team, led by
+ * David van der Spoel, Berk Hess, Erik Lindahl, and including many
+ * others, as listed in the AUTHORS file in the top-level source
+ * directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests for functionality of analysis data lifetime module.
+ *
+ * These tests check that gmx::AnalysisDataLifetimeModule computes lifetimes
+ * correctly with simple input data.
+ * Checking is done using gmx::test::AnalysisDataTestFixture and reference
+ * data. Also the input data is written to the reference data to catch
+ * out-of-date reference.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_analysisdata
+ */
+#include <gtest/gtest.h>
+
+#include "gromacs/analysisdata/analysisdata.h"
+#include "gromacs/analysisdata/modules/lifetime.h"
+
+#include "testutils/datatest.h"
+#include "testutils/testasserts.h"
+
+using gmx::test::AnalysisDataTestInput;
+
+namespace
+{
+
+// Simple input data for gmx::AnalysisDataLifetimeModule tests.
+class SimpleInputData
+{
+ public:
+ static const AnalysisDataTestInput &get()
+ {
+ static SimpleInputData singleton;
+ return singleton.data_;
+ }
+
+ SimpleInputData() : data_(1, false)
+ {
+ data_.setColumnCount(0, 3);
+ data_.addFrameWithValues(1.0, 1.0, 1.0, 1.0);
+ data_.addFrameWithValues(2.0, 1.0, 0.0, 1.0);
+ data_.addFrameWithValues(3.0, 0.0, 1.0, 1.0);
+ }
+
+ private:
+ AnalysisDataTestInput data_;
+};
+
+// Input data with multiple data sets for gmx::AnalysisDataLifetimeModule tests.
+class MultiDataSetInputData
+{
+ public:
+ static const AnalysisDataTestInput &get()
+ {
+ static MultiDataSetInputData singleton;
+ return singleton.data_;
+ }
+
+ MultiDataSetInputData() : data_(2, false)
+ {
+ using gmx::test::AnalysisDataTestInputFrame;
+ data_.setColumnCount(0, 2);
+ data_.setColumnCount(1, 2);
+ AnalysisDataTestInputFrame &frame1 = data_.addFrame(1.0);
+ frame1.addPointSetWithValues(0, 0, 1.0, 1.0);
+ frame1.addPointSetWithValues(1, 0, 0.0, 0.0);
+ AnalysisDataTestInputFrame &frame2 = data_.addFrame(2.0);
+ frame2.addPointSetWithValues(0, 0, 1.0, 0.0);
+ frame2.addPointSetWithValues(1, 0, 1.0, 0.0);
+ AnalysisDataTestInputFrame &frame3 = data_.addFrame(3.0);
+ frame3.addPointSetWithValues(0, 0, 1.0, 0.0);
+ frame3.addPointSetWithValues(1, 0, 1.0, 1.0);
+ }
+
+ private:
+ AnalysisDataTestInput data_;
+};
+
+
+/********************************************************************
+ * Tests for gmx::AnalysisDataLifetimeModule.
+ */
+
+//! Test fixture for gmx::AnalysisDataLifetimeModule.
+typedef gmx::test::AnalysisDataTestFixture LifetimeModuleTest;
+
+TEST_F(LifetimeModuleTest, BasicTest)
+{
+ const AnalysisDataTestInput &input = SimpleInputData::get();
+ gmx::AnalysisData data;
+ ASSERT_NO_THROW_GMX(setupDataObject(input, &data));
+
+ gmx::AnalysisDataLifetimeModulePointer module(
+ new gmx::AnalysisDataLifetimeModule);
+ module->setCumulative(false);
+ data.addModule(module);
+
+ ASSERT_NO_THROW_GMX(addStaticCheckerModule(input, &data));
+ ASSERT_NO_THROW_GMX(addReferenceCheckerModule("InputData", &data));
+ ASSERT_NO_THROW_GMX(addReferenceCheckerModule("Lifetime", module.get()));
+ ASSERT_NO_THROW_GMX(presentAllData(input, &data));
+}
+
+TEST_F(LifetimeModuleTest, CumulativeTest)
+{
+ const AnalysisDataTestInput &input = SimpleInputData::get();
+ gmx::AnalysisData data;
+ ASSERT_NO_THROW_GMX(setupDataObject(input, &data));
+
+ gmx::AnalysisDataLifetimeModulePointer module(
+ new gmx::AnalysisDataLifetimeModule);
+ module->setCumulative(true);
+ data.addModule(module);
+
+ ASSERT_NO_THROW_GMX(addStaticCheckerModule(input, &data));
+ ASSERT_NO_THROW_GMX(addReferenceCheckerModule("InputData", &data));
+ ASSERT_NO_THROW_GMX(addReferenceCheckerModule("Lifetime", module.get()));
+ ASSERT_NO_THROW_GMX(presentAllData(input, &data));
+}
+
+TEST_F(LifetimeModuleTest, HandlesMultipleDataSets)
+{
+ const AnalysisDataTestInput &input = MultiDataSetInputData::get();
+ gmx::AnalysisData data;
+ ASSERT_NO_THROW_GMX(setupDataObject(input, &data));
+
+ gmx::AnalysisDataLifetimeModulePointer module(
+ new gmx::AnalysisDataLifetimeModule);
+ module->setCumulative(false);
+ data.addModule(module);
+
+ ASSERT_NO_THROW_GMX(addStaticCheckerModule(input, &data));
+ ASSERT_NO_THROW_GMX(addReferenceCheckerModule("InputData", &data));
+ ASSERT_NO_THROW_GMX(addReferenceCheckerModule("Lifetime", module.get()));
+ ASSERT_NO_THROW_GMX(presentAllData(input, &data));
+}
+
+} // namespace
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <AnalysisData Name="InputData">
+ <DataFrame Name="Frame0">
+ <Real Name="X">1.000000</Real>
+ <DataValues>
+ <Int Name="Count">3</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">2.000000</Real>
+ <DataValues>
+ <Int Name="Count">3</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame2">
+ <Real Name="X">3.000000</Real>
+ <DataValues>
+ <Int Name="Count">3</Int>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
+ <AnalysisData Name="Lifetime">
+ <DataFrame Name="Frame0">
+ <Real Name="X">0.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">0.666667</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">1.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">0.500000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame2">
+ <Real Name="X">2.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <AnalysisData Name="InputData">
+ <DataFrame Name="Frame0">
+ <Real Name="X">1.000000</Real>
+ <DataValues>
+ <Int Name="Count">3</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">2.000000</Real>
+ <DataValues>
+ <Int Name="Count">3</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame2">
+ <Real Name="X">3.000000</Real>
+ <DataValues>
+ <Int Name="Count">3</Int>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
+ <AnalysisData Name="Lifetime">
+ <DataFrame Name="Frame0">
+ <Real Name="X">0.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">2.333333</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">1.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">1.500000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame2">
+ <Real Name="X">2.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <AnalysisData Name="InputData">
+ <DataFrame Name="Frame0">
+ <Real Name="X">1.000000</Real>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <Int Name="DataSet">0</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <Int Name="DataSet">1</Int>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">2.000000</Real>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <Int Name="DataSet">0</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ </DataValues>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <Int Name="DataSet">1</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame2">
+ <Real Name="X">3.000000</Real>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <Int Name="DataSet">0</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ </DataValues>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <Int Name="DataSet">1</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
+ <AnalysisData Name="Lifetime">
+ <DataFrame Name="Frame0">
+ <Real Name="X">0.000000</Real>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <DataValue>
+ <Real Name="Value">0.333333</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.333333</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">1.000000</Real>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.500000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame2">
+ <Real Name="X">2.000000</Real>
+ <DataValues>
+ <Int Name="Count">2</Int>
+ <DataValue>
+ <Real Name="Value">1.000000</Real>
+ </DataValue>
+ <DataValue>
+ <Real Name="Value">0.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
+</ReferenceData>
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013, by the GROMACS development team, led by
* David van der Spoel, Berk Hess, Erik Lindahl, and including many
* others, as listed in the AUTHORS file in the top-level source
* directory and at http://www.gromacs.org.
void TrajectoryAnalysisModule::registerBasicDataset(AbstractAnalysisData *data,
const char *name)
{
+ GMX_RELEASE_ASSERT(data != NULL, "Attempting to register NULL data");
// TODO: Strong exception safety should be possible to implement.
GMX_RELEASE_ASSERT(impl_->datasets_.find(name) == impl_->datasets_.end(),
"Duplicate data set name registered");
: TrajectoryAnalysisModule(name, shortDescription),
selOpt_(NULL),
bDump_(false), bTotNorm_(false), bFracNorm_(false), bResInd_(false),
- top_(NULL), occupancyModule_(new AnalysisDataAverageModule())
+ bCumulativeLifetimes_(true), top_(NULL),
+ occupancyModule_(new AnalysisDataAverageModule()),
+ lifetimeModule_(new AnalysisDataLifetimeModule())
{
registerAnalysisDataset(&sdata_, "size");
registerAnalysisDataset(&cdata_, "cfrac");
registerAnalysisDataset(&mdata_, "mask");
occupancyModule_->setXAxis(1.0, 1.0);
registerBasicDataset(occupancyModule_.get(), "occupancy");
+ registerBasicDataset(lifetimeModule_.get(), "lifetime");
}
"present, with [TT]maxsel[tt] all atoms possibly selected by the",
"selection are present, and with [TT]selected[tt] only atoms that are",
"selected at least in one frame are present.[PAR]",
- "With [TT]-om[tt], [TT]-of[tt] and [TT]-ofpdb[tt], only one selection",
- "can be provided. [TT]-om[tt] and [TT]-of[tt] only accept dynamic",
+ "With [TT]-olt[tt], a histogram is produced that shows the number of",
+ "selected positions as a function of the time the position was",
+ "continuously selected. [TT]-cumlt[tt] can be used to control whether",
+ "subintervals of longer intervals are included in the histogram.[PAR]",
+ "With [TT]-om[tt], [TT]-of[tt], [TT]-olt[tt], and [TT]-ofpdb[tt],",
+ "only one selection can be provided.",
+ "[TT]-om[tt], [TT]-of[tt], and [TT]-olt[tt] only accept dynamic",
"selections."
};
options->addOption(FileNameOption("ofpdb").filetype(eftPDB).outputFile()
.store(&fnPDB_).defaultBasename("occupancy")
.description("PDB file with occupied fraction for selected positions"));
+ options->addOption(FileNameOption("olt").filetype(eftPlot).outputFile()
+ .store(&fnLifetime_).defaultBasename("lifetime")
+ .description("Lifetime histogram"));
selOpt_ = options->addOption(SelectionOption("select").storeVector(&sel_)
.required().multiValue()
options->addOption(StringOption("pdbatoms").store(&pdbAtoms_)
.enumValue(cPDBAtomsEnum).defaultEnumIndex(0)
.description("Atoms to write with -ofpdb"));
+ options->addOption(BooleanOption("cumlt").store(&bCumulativeLifetimes_)
+ .description("Cumulate subintervals of longer intervals in -olt"));
}
void
settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
}
if ((!fnIndex_.empty() && bDump_)
- || !fnMask_.empty() || !fnOccupancy_.empty() || !fnPDB_.empty())
+ || !fnMask_.empty() || !fnOccupancy_.empty() || !fnPDB_.empty()
+ || !fnLifetime_.empty())
{
selOpt_->setValueCount(1);
}
Select::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
{
- if (!sel_[0].isDynamic() && (!fnMask_.empty() || !fnOccupancy_.empty()))
+ if (!sel_[0].isDynamic()
+ && (!fnMask_.empty() || !fnOccupancy_.empty() || !fnLifetime_.empty()))
{
GMX_THROW(InconsistentInputError(
- "-om or -of are not meaningful with a static selection"));
+ "-om, -of, or -olt are not meaningful with a static selection"));
}
bResInd_ = (resNumberType_ == "index");
idata_.addModule(writer);
}
+ // TODO: Write out the mask and lifetimes for all the selections.
mdata_.setColumnCount(0, sel_[0].posCount());
mdata_.addModule(occupancyModule_);
+ mdata_.addModule(lifetimeModule_);
+ lifetimeModule_->setCumulative(bCumulativeLifetimes_);
if (!fnMask_.empty())
{
AnalysisDataPlotModulePointer plot(
plot->setYLabel("Occupied fraction");
occupancyModule_->addColumnModule(0, 1, plot);
}
+ if (!fnLifetime_.empty())
+ {
+ AnalysisDataPlotModulePointer plot(
+ new AnalysisDataPlotModule(settings.plotSettings()));
+ plot->setFileName(fnLifetime_);
+ plot->setTitle("Lifetime histogram");
+ plot->setXAxisIsTime();
+ plot->setYLabel("Number of occurrences");
+ lifetimeModule_->addModule(plot);
+ }
top_ = ⊤
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2011,2012, by the GROMACS development team, led by
+ * Copyright (c) 2011,2012,2013, by the GROMACS development team, led by
* David van der Spoel, Berk Hess, Erik Lindahl, and including many
* others, as listed in the AUTHORS file in the top-level source
* directory and at http://www.gromacs.org.
#include "../analysismodule.h"
#include "gromacs/analysisdata/analysisdata.h"
#include "gromacs/analysisdata/modules/average.h"
+#include "gromacs/analysisdata/modules/lifetime.h"
#include "gromacs/selection/selection.h"
namespace gmx
virtual void writeOutput();
private:
- SelectionList sel_;
- SelectionOptionInfo *selOpt_;
-
- std::string fnSize_;
- std::string fnFrac_;
- std::string fnIndex_;
- std::string fnNdx_;
- std::string fnMask_;
- std::string fnOccupancy_;
- std::string fnPDB_;
- bool bDump_;
- bool bTotNorm_;
- bool bFracNorm_;
- bool bResInd_;
- std::string resNumberType_;
- std::string pdbAtoms_;
-
- const TopologyInformation *top_;
- std::vector<int> totsize_;
- AnalysisData sdata_;
- AnalysisData cdata_;
- AnalysisData idata_;
- AnalysisData mdata_;
- AnalysisDataAverageModulePointer occupancyModule_;
+ SelectionList sel_;
+ SelectionOptionInfo *selOpt_;
+
+ std::string fnSize_;
+ std::string fnFrac_;
+ std::string fnIndex_;
+ std::string fnNdx_;
+ std::string fnMask_;
+ std::string fnOccupancy_;
+ std::string fnPDB_;
+ std::string fnLifetime_;
+ bool bDump_;
+ bool bTotNorm_;
+ bool bFracNorm_;
+ bool bResInd_;
+ bool bCumulativeLifetimes_;
+ std::string resNumberType_;
+ std::string pdbAtoms_;
+
+ const TopologyInformation *top_;
+ std::vector<int> totsize_;
+ AnalysisData sdata_;
+ AnalysisData cdata_;
+ AnalysisData idata_;
+ AnalysisData mdata_;
+ AnalysisDataAverageModulePointer occupancyModule_;
+ AnalysisDataLifetimeModulePointer lifetimeModule_;
};
} // namespace analysismodules
</DataValues>
</DataFrame>
</AnalysisData>
+ <AnalysisData Name="lifetime">
+ <DataFrame Name="Frame0">
+ <Real Name="X">0.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">8.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ <DataFrame Name="Frame1">
+ <Real Name="X">0.000000</Real>
+ <DataValues>
+ <Int Name="Count">1</Int>
+ <DataValue>
+ <Real Name="Value">3.000000</Real>
+ </DataValue>
+ </DataValues>
+ </DataFrame>
+ </AnalysisData>
<AnalysisData Name="mask">
<DataFrame Name="Frame0">
<Real Name="X">0.000000</Real>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff