mrcdensitymapheader.cpp
readinp.cpp
fileioxdrserializer.cpp
+ xvgio.cpp
)
if (GMX_USE_TNG)
list(APPEND test_sources tngio.cpp)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and 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 tests for xvg file operations
+ *
+ * \author Joe Jordan <ejjordan@kth.se>
+ */
+#include "gmxpre.h"
+
+#include <numeric>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/textwriter.h"
+
+#include "testutils/testfilemanager.h"
+#include "testutils/testoptions.h"
+
+namespace gmx
+{
+namespace test
+{
+
+class XvgioTester
+{
+public:
+ XvgioTester() { referenceFilename_ = fileManager_.getTemporaryFilePath("ref.xvg"); }
+
+ const std::string& referenceFilename() const { return referenceFilename_; }
+
+ const std::string& referenceContents() const { return referenceContents_; }
+
+ void useStringAsXvgFile(const std::string& xvgString) { referenceContents_ = xvgString; }
+
+ void writeXvgFile()
+ {
+ gmx::TextWriter::writeFileFromString(referenceFilename(), referenceContents());
+ }
+
+ void compareValues(basic_mdspan<const double, dynamicExtents2D> ref,
+ basic_mdspan<const double, dynamicExtents2D> test)
+ {
+ // The xvg reading routines use a column-major layout, while we would
+ // like to enforce row major behaviour everywhere else. This requires
+ // this test to swap the orders between reference data and test data.
+ // Hence, we compare extent(0) with extent(1) and [i][j] with [j][i].
+ EXPECT_EQ(ref.extent(0), test.extent(1));
+ EXPECT_EQ(ref.extent(1), test.extent(0));
+
+ for (std::ptrdiff_t i = 0; i < ref.extent(0); i++)
+ {
+ for (std::ptrdiff_t j = 0; j < ref.extent(1); j++)
+ {
+ EXPECT_DOUBLE_EQ(ref[i][j], test[j][i]);
+ }
+ }
+ }
+
+private:
+ gmx::test::TestFileManager fileManager_;
+ std::string referenceFilename_;
+ std::string referenceContents_;
+};
+
+TEST(XvgioTest, readXvgIntWorks)
+{
+ XvgioTester xvgioTester;
+ xvgioTester.useStringAsXvgFile(
+ "1 2 3\n"
+ "4 5 6\n");
+ xvgioTester.writeXvgFile();
+ MultiDimArray<std::vector<double>, dynamicExtents2D> xvgTestData =
+ readXvgData(xvgioTester.referenceFilename());
+
+ const int numRows = 2;
+ const int numColumns = 3;
+ MultiDimArray<std::vector<double>, dynamicExtents2D> xvgRefData(numRows, numColumns);
+ std::iota(begin(xvgRefData), end(xvgRefData), 1);
+
+ xvgioTester.compareValues(xvgRefData.asConstView(), xvgTestData.asConstView());
+}
+
+TEST(XvgioTest, readXvgRealWorks)
+{
+ XvgioTester xvgioTester;
+ xvgioTester.useStringAsXvgFile(
+ "1.1 2.2\n"
+ "3.3 4.4\n"
+ "5.5 6.6\n");
+ xvgioTester.writeXvgFile();
+ MultiDimArray<std::vector<double>, dynamicExtents2D> xvgTestData =
+ readXvgData(xvgioTester.referenceFilename());
+
+ const int numRows = 3;
+ const int numColumns = 2;
+ MultiDimArray<std::vector<double>, dynamicExtents2D> xvgRefData(numRows, numColumns);
+ std::generate(begin(xvgRefData), end(xvgRefData), [n = 0.0]() mutable {
+ n += 1.1;
+ return n;
+ });
+ xvgioTester.compareValues(xvgRefData.asConstView(), xvgTestData.asConstView());
+}
+
+TEST(XvgioTest, readXvgIgnoreCommentLineWorks)
+{
+ XvgioTester xvgioTester;
+ xvgioTester.useStringAsXvgFile(
+ "1 2 3\n"
+ "#comment\n"
+ "4 5 6\n");
+ xvgioTester.writeXvgFile();
+
+ MultiDimArray<std::vector<double>, dynamicExtents2D> xvgTestData =
+ readXvgData(xvgioTester.referenceFilename());
+
+ const int numRows = 2;
+ const int numColumns = 3;
+ MultiDimArray<std::vector<double>, dynamicExtents2D> xvgRefData(numRows, numColumns);
+ std::iota(begin(xvgRefData), end(xvgRefData), 1);
+
+ xvgioTester.compareValues(xvgRefData.asConstView(), xvgTestData.asConstView());
+}
+
+// Todo: Remove this test once all calls to read_xvg have been ported to readXvgData
+TEST(XvgioTest, readXvgDeprecatedWorks)
+{
+ XvgioTester xvgioTester;
+ xvgioTester.useStringAsXvgFile(
+ "1 2 3\n"
+ "4 5 6\n");
+ xvgioTester.writeXvgFile();
+ std::vector<std::vector<double>> xvgData = { { 1, 4 }, { 2, 5 }, { 3, 6 } };
+
+ double** xvgTestData = nullptr;
+ int testNumColumns;
+ int testNumRows = read_xvg(xvgioTester.referenceFilename().c_str(), &xvgTestData, &testNumColumns);
+
+ double** xvgRefData = nullptr;
+ int refNumColumns = 3;
+ int refNumRows = 2;
+
+ EXPECT_EQ(refNumColumns, testNumColumns);
+ EXPECT_EQ(refNumRows, testNumRows);
+
+ // Set the reference data
+ snew(xvgRefData, refNumColumns);
+ for (int column = 0; column < refNumColumns; column++)
+ {
+ snew(xvgRefData[column], refNumRows);
+ for (int row = 0; row < refNumRows; row++)
+ {
+ xvgRefData[column][row] = xvgData[column][row];
+ }
+ }
+
+ // Check that the reference and test data match
+ for (int column = 0; column < refNumColumns; column++)
+ {
+ for (int row = 0; row < refNumRows; row++)
+ {
+ EXPECT_EQ(xvgRefData[column][row], xvgTestData[column][row]);
+ }
+ }
+
+ // Free the reference and test data memory
+ for (int column = 0; column < refNumColumns; column++)
+ {
+ sfree(xvgRefData[column]);
+ sfree(xvgTestData[column]);
+ }
+ sfree(xvgRefData);
+ sfree(xvgTestData);
+}
+
+} // namespace test
+} // namespace gmx
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/oenv.h"
#include "gromacs/math/vec.h"
+#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/binaryinformation.h"
#include "gromacs/utility/coolstuff.h"
#include "gromacs/utility/cstringutil.h"
int read_xvg(const char* fn, double*** y, int* ny)
{
- return read_xvg_legend(fn, y, ny, nullptr, nullptr);
+ gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgData =
+ readXvgData(std::string(fn));
+
+ int numColumns = xvgData.extent(0);
+ int numRows = xvgData.extent(1);
+
+ double** yy = nullptr;
+ snew(yy, numColumns);
+ for (int column = 0; column < numColumns; column++)
+ {
+ snew(yy[column], numRows);
+ for (int row = 0; row < numRows; row++)
+ {
+ yy[column][row] = xvgData.asConstView()[column][row];
+ }
+ }
+
+ *y = yy;
+ *ny = numColumns;
+ int nx = numRows;
+ return nx;
+}
+
+gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> readXvgData(const std::string& fn)
+{
+ FILE* fp = gmx_fio_fopen(fn.c_str(), "r");
+ char* ptr;
+ char* base = nullptr;
+ char* fmt = nullptr;
+ char* tmpbuf;
+ int len = STRLEN;
+
+ //! This only gets properly set after the first line of data is read
+ int numColumns = 0;
+ int numRows = 0;
+ snew(tmpbuf, len);
+ std::vector<double> xvgData;
+
+ for (int line = 0; (ptr = fgets3(fp, &tmpbuf, &len, 10 * STRLEN)) != nullptr && ptr[0] != '&'; ++line)
+ {
+ trim(ptr);
+ if (ptr[0] == '@' || ptr[0] == '#')
+ {
+ continue;
+ }
+ ++numRows;
+ if (numColumns == 0)
+ {
+ numColumns = wordcount(ptr);
+ if (numColumns == 0)
+ {
+ return {}; // There are no columns and hence no data to process
+ }
+ snew(fmt, 3 * numColumns + 1);
+ snew(base, 3 * numColumns + 1);
+ }
+ /* Initiate format string */
+ fmt[0] = '\0';
+ base[0] = '\0';
+ int columnCount = 0;
+ for (columnCount = 0; (columnCount < numColumns); columnCount++)
+ {
+ double lf;
+ std::strcpy(fmt, base);
+ std::strcat(fmt, "%lf");
+ int rval = sscanf(ptr, fmt, &lf);
+ if ((rval == EOF) || (rval == 0))
+ {
+ break;
+ }
+ xvgData.push_back(lf);
+ srenew(fmt, 3 * (numColumns + 1) + 1);
+ srenew(base, 3 * numColumns + 1);
+ std::strcat(base, "%*s");
+ }
+
+ if (columnCount != numColumns)
+ {
+ fprintf(stderr, "Only %d columns on line %d in file %s\n", columnCount, line, fn.c_str());
+ for (; (columnCount < numColumns); columnCount++)
+ {
+ xvgData.push_back(0.0);
+ }
+ }
+ }
+ gmx_fio_fclose(fp);
+
+ sfree(tmpbuf);
+ sfree(base);
+ sfree(fmt);
+
+ gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgDataAsArray(numRows, numColumns);
+ std::copy(std::begin(xvgData), std::end(xvgData), begin(xvgDataAsArray.asView()));
+
+ gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgDataAsArrayTransposed(
+ numColumns, numRows);
+ for (std::ptrdiff_t row = 0; row < numRows; ++row)
+ {
+ for (std::ptrdiff_t column = 0; column < numColumns; ++column)
+ {
+ xvgDataAsArrayTransposed(column, row) = xvgDataAsArray(row, column);
+ }
+ }
+
+ return xvgDataAsArrayTransposed;
}
void write_xvg(const char* fn, const char* title, int nx, int ny, real** y, const char** leg, const gmx_output_env_t* oenv)
#include <string>
#include <vector>
+#include "gromacs/math/multidimarray.h"
+#include "gromacs/mdspan/extensions.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
* 0 is the first y legend, the legend string will be NULL when not present.
*/
+/* \brief Read only the data from an xvg file for post processing.
+ *
+ * Note: this function is deprecated in favor of readXvg, which is
+ * used under the hood in this function.
+ *
+ * \param[out] nx Number of rows.
+ * \param[in] fn Xvg file to read.
+ * \param[in/out] y Pointer to 2D array (allocated by the routine).
+ * \param[in/out] ny Number of columns.
+ *
+ * Todo: Port all read_xvg calls to use readXvgData
+ */
int read_xvg(const char* fn, double*** y, int* ny);
-/* As read_xvg_legend, but does not read legends. */
+
+/* \brief Read only the data from an xvg file for post processing.
+ *
+ * \param[out] XvgData Data in row major.
+ * \param[in] fn Xvg file to read.
+ */
+gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> readXvgData(const std::string& fn);
+
void write_xvg(const char* fn,
const char* title,
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
return basicMdspan.data() + basicMdspan.mapping().required_span_size();
}
+//! Convenience type for often-used two dimensional extents
+using dynamicExtents2D = extents<dynamic_extent, dynamic_extent>;
+
//! Convenience type for often-used three dimensional extents
using dynamicExtents3D = extents<dynamic_extent, dynamic_extent, dynamic_extent>;
#include "gromacs/fileio/xvgr.h"
#include "gromacs/math/functions.h"
+#include "gromacs/math/multidimarray.h"
#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
+#include "gromacs/mdspan/extensions.h"
#include "gromacs/mdtypes/fcdata.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/nblist.h"
static void read_tables(FILE* fp, const char* filename, int ntab, int angle, t_tabledata td[])
{
char buf[STRLEN];
- double **yy = nullptr, start, end, dx0, dx1, ssd, vm, vp, f, numf;
- int k, i, nx, nx0 = 0, ny, nny, ns;
+ double start, end, dx0, dx1, ssd, vm, vp, f, numf;
+ int k, i, nx0 = 0, nny, ns;
gmx_bool bAllZero, bZeroV, bZeroF;
double tabscale;
nny = 2 * ntab + 1;
std::string libfn = gmx::findLibraryFile(filename);
- nx = read_xvg(libfn.c_str(), &yy, &ny);
- if (ny != nny)
+ gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgData = readXvgData(libfn);
+ int numColumns = xvgData.extent(0);
+ if (numColumns != nny)
{
gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d", libfn.c_str(),
- ny, nny);
+ numColumns, nny);
}
+ int numRows = xvgData.extent(1);
+
+ const auto& yy = xvgData.asView();
if (angle == 0)
{
if (yy[0][0] != 0.0)
start = -180.0;
}
end = 180.0;
- if (yy[0][0] != start || yy[0][nx - 1] != end)
+ if (yy[0][0] != start || yy[0][numRows - 1] != end)
{
gmx_fatal(FARGS, "The angles in file %s should go from %f to %f instead of %f to %f\n",
- libfn.c_str(), start, end, yy[0][0], yy[0][nx - 1]);
+ libfn.c_str(), start, end, yy[0][0], yy[0][numRows - 1]);
}
}
- tabscale = (nx - 1) / (yy[0][nx - 1] - yy[0][0]);
+ tabscale = (numRows - 1) / (yy[0][numRows - 1] - yy[0][0]);
if (fp)
{
- fprintf(fp, "Read user tables from %s with %d data points.\n", libfn.c_str(), nx);
+ fprintf(fp, "Read user tables from %s with %d data points.\n", libfn.c_str(), numRows);
if (angle == 0)
{
fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
{
bZeroV = TRUE;
bZeroF = TRUE;
- for (i = 0; (i < nx); i++)
+ for (i = 0; (i < numRows); i++)
{
if (i >= 2)
{
if (!bZeroV && bZeroF)
{
- set_forces(fp, angle, nx, 1 / tabscale, yy[1 + k * 2], yy[1 + k * 2 + 1], k);
+ set_forces(fp, angle, numRows, 1 / tabscale, yy[1 + k * 2].data(),
+ yy[1 + k * 2 + 1].data(), k);
}
else
{
*/
ssd = 0;
ns = 0;
- for (i = 1; (i < nx - 1); i++)
+ for (i = 1; (i < numRows - 1); i++)
{
vm = yy[1 + 2 * k][i - 1];
vp = yy[1 + 2 * k][i + 1];
for (k = 0; (k < ntab); k++)
{
- init_table(nx, nx0, tabscale, &(td[k]), TRUE);
- for (i = 0; (i < nx); i++)
+ init_table(numRows, nx0, tabscale, &(td[k]), TRUE);
+ for (i = 0; (i < numRows); i++)
{
td[k].x[i] = yy[0][i];
td[k].v[i] = yy[2 * k + 1][i];
td[k].f[i] = yy[2 * k + 2][i];
}
}
- for (i = 0; (i < ny); i++)
- {
- sfree(yy[i]);
- }
- sfree(yy);
}
static void done_tabledata(t_tabledata* td)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff