};
/*! \brief Routine to compute ACF using FFT. */
-static void low_do_four_core(int nfour, int nframes, real c1[], real cfour[],
+static void low_do_four_core(int nframes, real c1[], real cfour[],
int nCos)
{
int i = 0;
-
+ std::vector<std::vector<real> > data;
+ data.resize(1);
+ data[0].resize(nframes, 0);
switch (nCos)
{
case enNorm:
for (i = 0; (i < nframes); i++)
{
- cfour[i] = c1[i];
+ data[0][i] = c1[i];
}
break;
case enCos:
for (i = 0; (i < nframes); i++)
{
- cfour[i] = cos(c1[i]);
+ data[0][i] = cos(c1[i]);
}
break;
case enSin:
for (i = 0; (i < nframes); i++)
{
- cfour[i] = sin(c1[i]);
+ data[0][i] = sin(c1[i]);
}
break;
default:
gmx_fatal(FARGS, "nCos = %d, %s %d", nCos, __FILE__, __LINE__);
}
- many_auto_correl(1, nframes, nfour, &cfour);
+ many_auto_correl(&data);
+ for (i = 0; (i < nframes); i++)
+ {
+ cfour[i] = data[0][i];
+ }
}
/*! \brief Routine to comput ACF without FFT. */
}
/*! \brief High level ACF routine. */
-void do_four_core(unsigned long mode, int nfour, int nf2, int nframes,
- real c1[], real csum[], real ctmp[])
+static void do_four_core(unsigned long mode, int nframes,
+ real c1[], real csum[], real ctmp[])
{
real *cfour;
char buf[32];
real fac;
int j, m, m1;
- snew(cfour, nfour);
+ snew(cfour, nframes);
if (MODE(eacNormal))
{
/********************************************
* N O R M A L
********************************************/
- low_do_four_core(nfour, nf2, c1, csum, enNorm);
+ low_do_four_core(nframes, c1, csum, enNorm);
}
else if (MODE(eacCos))
{
/* Copy the data to temp array. Since we need it twice
* we can't overwrite original.
*/
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
ctmp[j] = c1[j];
}
/* Cosine term of AC function */
- low_do_four_core(nfour, nf2, ctmp, cfour, enCos);
- for (j = 0; (j < nf2); j++)
+ low_do_four_core(nframes, ctmp, cfour, enCos);
+ for (j = 0; (j < nframes); j++)
{
c1[j] = cfour[j];
}
/* Sine term of AC function */
- low_do_four_core(nfour, nf2, ctmp, cfour, enSin);
- for (j = 0; (j < nf2); j++)
+ low_do_four_core(nframes, ctmp, cfour, enSin);
+ for (j = 0; (j < nframes); j++)
{
c1[j] += cfour[j];
csum[j] = c1[j];
/* Because of normalization the number of -0.5 to subtract
* depends on the number of data points!
*/
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
- csum[j] = -0.5*(nf2-j);
+ csum[j] = -0.5*(nframes-j);
}
/***** DIAGONAL ELEMENTS ************/
for (m = 0; (m < DIM); m++)
{
/* Copy the vector data in a linear array */
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
ctmp[j] = gmx::square(c1[DIM*j+m]);
}
if (debug)
{
sprintf(buf, "c1diag%d.xvg", m);
- dump_tmp(buf, nf2, ctmp);
+ dump_tmp(buf, nframes, ctmp);
}
- low_do_four_core(nfour, nf2, ctmp, cfour, enNorm);
+ low_do_four_core(nframes, ctmp, cfour, enNorm);
if (debug)
{
sprintf(buf, "c1dfout%d.xvg", m);
- dump_tmp(buf, nf2, cfour);
+ dump_tmp(buf, nframes, cfour);
}
fac = 1.5;
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
csum[j] += fac*(cfour[j]);
}
{
/* Copy the vector data in a linear array */
m1 = (m+1) % DIM;
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
ctmp[j] = c1[DIM*j+m]*c1[DIM*j+m1];
}
if (debug)
{
sprintf(buf, "c1off%d.xvg", m);
- dump_tmp(buf, nf2, ctmp);
+ dump_tmp(buf, nframes, ctmp);
}
- low_do_four_core(nfour, nf2, ctmp, cfour, enNorm);
+ low_do_four_core(nframes, ctmp, cfour, enNorm);
if (debug)
{
sprintf(buf, "c1ofout%d.xvg", m);
- dump_tmp(buf, nf2, cfour);
+ dump_tmp(buf, nframes, cfour);
}
fac = 3.0;
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
csum[j] += fac*cfour[j];
}
* First for XX, then for YY, then for ZZ
* After that we sum them and normalise
*/
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
csum[j] = 0.0;
}
for (m = 0; (m < DIM); m++)
{
/* Copy the vector data in a linear array */
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
ctmp[j] = c1[DIM*j+m];
}
- low_do_four_core(nfour, nf2, ctmp, cfour, enNorm);
- for (j = 0; (j < nf2); j++)
+ low_do_four_core(nframes, ctmp, cfour, enNorm);
+ for (j = 0; (j < nframes); j++)
{
csum[j] += cfour[j];
}
}
sfree(cfour);
- for (j = 0; (j < nf2); j++)
+ for (j = 0; (j < nframes); j++)
{
c1[j] = csum[j]/(real)(nframes-j);
}
int eFitFn)
{
FILE *fp, *gp = NULL;
- int i, nfour;
+ int i;
real *csum;
real *ctmp, *fit;
real sum, Ct2av, Ctav;
gmx::boolToString(bNormalize));
printf("mode = %lu, dt = %g, nrestart = %d\n", mode, dt, nrestart);
}
- if (bFour)
- {
- /* For FTT corr., we need to pad the data with at least nframes zeros */
- nfour = 2;
- while (2*nframes > nfour)
- {
- nfour *= 2;
- }
- if (debug)
- {
- fprintf(debug, "Using FFT to calculate %s, #points for FFT = %d\n",
- title, nfour);
- }
-
- /* Allocate temp arrays */
- snew(csum, nfour);
- snew(ctmp, nfour);
- }
- else
- {
- nfour = 0; /* To keep the compiler happy */
- snew(csum, nframes);
- snew(ctmp, nframes);
- }
+ /* Allocate temp arrays */
+ snew(csum, nframes);
+ snew(ctmp, nframes);
/* Loop over items (e.g. molecules or dihedrals)
* In this loop the actual correlation functions are computed, but without
if (bFour)
{
- do_four_core(mode, nfour, nframes, nframes, c1[i], csum, ctmp);
+ do_four_core(mode, nframes, c1[i], csum, ctmp);
}
else
{
#include "manyautocorrelation.h"
-#include <stdio.h>
-#include <stdlib.h>
-
-#include <cmath>
-
#include <algorithm>
#include "gromacs/fft/fft.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/gmxomp.h"
-#include "gromacs/utility/smalloc.h"
-int many_auto_correl(int nfunc, int ndata, int nfft, real **c)
+int many_auto_correl(std::vector<std::vector<real> > *c)
{
+ size_t nfunc = (*c).size();
+ if (nfunc == 0)
+ {
+ GMX_THROW(gmx::InconsistentInputError("Empty array of vectors supplied"));
+ }
+ size_t ndata = (*c)[0].size();
+ if (ndata == 0)
+ {
+ GMX_THROW(gmx::InconsistentInputError("Empty vector supplied"));
+ }
+#ifndef NDEBUG
+ for (size_t i = 1; i < nfunc; i++)
+ {
+ if ((*c)[i].size() != ndata)
+ {
+ char buf[256];
+ snprintf(buf, sizeof(buf), "Vectors of different lengths supplied (%d %d)",
+ static_cast<int>((*c)[i].size()),
+ static_cast<int>(ndata));
+ GMX_THROW(gmx::InconsistentInputError(buf));
+ }
+ }
+#endif
+ // Add buffer size to the arrays.
+ size_t nfft = (3*ndata/2) + 1;
+ // Pad arrays with zeros
+ for (auto &i : *c)
+ {
+ i.resize(nfft, 0);
+ }
#pragma omp parallel
{
try
{
- typedef real complex[2];
- int i, j;
- gmx_fft_t fft1;
- complex *in, *out;
- int i0, i1;
- int nthreads, thread_id;
+ gmx_fft_t fft1;
+ std::vector<real> in, out;
- nthreads = gmx_omp_get_max_threads();
- thread_id = gmx_omp_get_thread_num();
- if ((0 == thread_id))
- {
- // fprintf(stderr, "There are %d threads for correlation functions\n", nthreads);
- }
- i0 = thread_id*nfunc/nthreads;
- i1 = std::min(nfunc, (thread_id+1)*nfunc/nthreads);
+ int nthreads = gmx_omp_get_max_threads();
+ int thread_id = gmx_omp_get_thread_num();
+ int i0 = (thread_id*nfunc)/nthreads;
+ int i1 = std::min(nfunc, ((thread_id+1)*nfunc)/nthreads);
gmx_fft_init_1d(&fft1, nfft, GMX_FFT_FLAG_CONSERVATIVE);
/* Allocate temporary arrays */
- snew(in, nfft);
- snew(out, nfft);
- for (i = i0; (i < i1); i++)
+ in.resize(2*nfft, 0);
+ out.resize(2*nfft, 0);
+ for (int i = i0; (i < i1); i++)
{
- for (j = 0; j < ndata; j++)
- {
- in[j][0] = c[i][j];
- in[j][1] = 0;
- }
- for (; (j < nfft); j++)
+ for (size_t j = 0; j < ndata; j++)
{
- in[j][0] = in[j][1] = 0;
+ in[2*j+0] = (*c)[i][j];
+ in[2*j+1] = 0;
}
-
- gmx_fft_1d(fft1, GMX_FFT_BACKWARD, (void *)in, (void *)out);
- for (j = 0; j < nfft; j++)
+ gmx_fft_1d(fft1, GMX_FFT_BACKWARD, (void *)in.data(), (void *)out.data());
+ for (size_t j = 0; j < nfft; j++)
{
- in[j][0] = (out[j][0]*out[j][0] + out[j][1]*out[j][1])/nfft;
- in[j][1] = 0;
+ in[2*j+0] = (out[2*j+0]*out[2*j+0] + out[2*j+1]*out[2*j+1])/nfft;
+ in[2*j+1] = 0;
}
- for (; (j < nfft); j++)
+ gmx_fft_1d(fft1, GMX_FFT_FORWARD, (void *)in.data(), (void *)out.data());
+ for (size_t j = 0; (j < nfft); j++)
{
- in[j][0] = in[j][1] = 0;
- }
-
- gmx_fft_1d(fft1, GMX_FFT_FORWARD, (void *)in, (void *)out);
- for (j = 0; (j < nfft); j++)
- {
- c[i][j] = out[j][0];
+ (*c)[i][j] = out[2*j+0];
}
}
/* Free the memory */
gmx_fft_destroy(fft1);
- sfree(in);
- sfree(out);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
}
- // gmx_fft_cleanup();
+ for (auto &i : *c)
+ {
+ i.resize(ndata);
+ }
+
return 0;
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014, by the GROMACS development team, led by
+ * Copyright (c) 2014,2016, 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.
#ifndef GMX_MANYAUTOCORRELATION_H
#define GMX_MANYAUTOCORRELATION_H
+#include <vector>
+
#include "gromacs/fft/fft.h"
#include "gromacs/utility/real.h"
-#ifdef __cplusplus
-extern "C" {
-#endif
-
/*! \brief
* Perform many autocorrelation calculations.
*
* This routine performs many autocorrelation function calculations using FFTs.
- * The GROMACS FFT library wrapper is employed. On return the c[] arrays contain
+ * The GROMACS FFT library wrapper is employed. On return the c vector contain
* a symmetric function that is useful for further FFT:ing, for instance in order to
* compute spectra.
*
+ * The vectors c[i] should all have the same length, but this is not checked for.
+ *
+ * The c arrays will be extend and filled with zero beyond ndata before
+ * computing the correlation.
+ *
* The functions uses OpenMP parallellization.
*
- * \param[in] nfunc Number of data functions to autocorrelate
- * \param[in] ndata Number of valid data points in the data
- * \param[in] nfft Length of the data arrays, this should at least be 50% larger than ndata. The c arrays will filled with zero beyond ndata before computing the correlation.
- * \param[inout] c Data array of size nfunc x nfft, will also be used for output
+ * \param[inout] c Data array
* \return fft error code, or zero if everything went fine (see fft/fft.h)
+ * \throws gmx::InconsistentInputError if the input is inconsistent.
*/
-int many_auto_correl(int nfunc, int ndata, int nfft, real **c);
-
-#ifdef __cplusplus
-}
-#endif
+int many_auto_correl(std::vector<std::vector<real> > *c);
#endif
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2014, by the GROMACS development team, led by
+# Copyright (c) 2014,2016, 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.
gmx_add_unit_test(CorrelationsTest correlations-test
autocorr.cpp
+ manyautocorrelation.cpp
correlationdataset.cpp
expfit.cpp)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014,2015,2016, 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 low level test of manyautocorrelation routines
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ * \ingroup module_correlationfunctions
+ */
+#include "gmxpre.h"
+
+#include "gromacs/correlationfunctions/manyautocorrelation.h"
+
+#include <cmath>
+
+#include <memory>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/utility/exceptions.h"
+//#include "gromacs/utility/gmxassert.h"
+
+#include "testutils/testasserts.h"
+#include "testutils/testfilemanager.h"
+
+namespace gmx
+{
+namespace
+{
+
+class ManyAutocorrelationTest : public ::testing::Test
+{
+};
+
+TEST_F (ManyAutocorrelationTest, Empty)
+{
+ std::vector<std::vector<real> > c;
+ EXPECT_THROW_GMX(many_auto_correl(&c), gmx::InconsistentInputError);
+}
+
+TEST_F (ManyAutocorrelationTest, DifferentLength)
+{
+ std::vector<std::vector<real> > c;
+ c.resize(3);
+ c[0].resize(10);
+ c[1].resize(10);
+ c[2].resize(8);
+ EXPECT_THROW_GMX(many_auto_correl(&c), gmx::InconsistentInputError);
+}
+
+}
+
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff