%
% This file is part of the GROMACS molecular simulation package.
%
-% Copyright (c) 2013,2014, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2015, 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.
actually the values do not converge.
%} % Brace matches ifthenelse test for gmxlite
+\section{Curve fitting in \gromacs}
+\subsection{Sum of exponential functions}
+Sometimes it is useful to fit a curve to an analytical function, for
+example in the case of autocorrelation functions with noisy
+tails. {\gromacs} is not a general purpose curve-fitting tool however
+and therefore {\gromacs} only supports a limited number of
+functions.
+Table~\ref{tab:fitfn} lists the available options with the
+corresponding command-line options. The underlying routines for
+fitting use the Levenberg-Marquardt algorithm as implemented in the
+{\tt lmfit} package~\cite{lmfit} (a bare-bones version of which is
+included in {\gromacs} in which an option for error-weighted fitting
+was implemented).
+\begin{table}[ht]
+\centering
+\caption{Overview of fitting functions supported in (most) analysis tools
+ that compute autocorrelation functions. The ``Note'' column describes
+ properties of the output parameters.}
+\label{tab:fitfn}
+\begin{tabular}{lcc}
+\hline
+Command & Functional form $f(t)$& Note\\
+line option & & \\
+\hline
+exp & $e^{-t/{a_0}}$ &\\
+aexp & $a_1e^{-t/{a_0}}$ &\\
+exp_exp & $a_1e^{-t/{a_0}}+(1-a_1)e^{-t/{a_2}}$ & $a_2\ge a_0\ge 0$\\
+exp5 & $a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_4$ &$a_2\ge a_0\ge 0$\\
+exp7 &
+$a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_6$&$a_4\ge a_2\ge
+a_0 \ge0$\\
+exp9 &
+$a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_7e^{-t/{a_6}}+a_8$&$a_6\ge
+a_4\ge a_2\ge a_0\ge 0$\\
+\hline
+\end{tabular}
+\end{table}
+
+\subsection{Error estimation}
+Under the hood {\gromacs} implements some more fitting functions,
+namely a function to estimate the error in time-correlated data due to Hess~\cite{Hess2002a}:
+\beq
+\varepsilon^2(t) =
+\alpha\tau_1\left(1+\frac{\tau_1}{t}\left(e^{-t/\tau_1}-1\right)\right)
+ + (1-\alpha)\tau_2\left(1+\frac{\tau_2}{t}\left(e^{-t/\tau_2}-1\right)\right)
+\eeq
+where $\tau_1$ and
+$\tau_2$ are time constants (with $\tau_2 \ge \tau_1$) and $\alpha$
+usually is close to 1 (in the fitting procedure it is enforced that
+$0\leq\alpha\leq 1$). This is used in {\tt gmx analyze} for error
+estimation using
+\beq
+\lim_{t\rightarrow\infty}\varepsilon(t) = \sigma\sqrt{\frac{2(\alpha\tau_1+(1-\alpha)\tau_2)}{T}}
+\eeq
+where $\sigma$ is the standard deviation of the data set and $T$ is
+the total simulation time~\cite{Hess2002a}.
+
+\subsection{Interphase boundary demarcation}
+In order to determine the position and width of an interface,
+Steen-S{\ae}thre {\em et al.} fitted a density profile to the
+following function
+\beq
+f(x) ~=~ \frac{a_0+a_1}{2} - \frac{a_0-a_1}{2}{\rm
+ erf}\left(\frac{x-a_2}{a_3^2}\right)
+\eeq
+where $a_0$ and $a_1$ are densities of different phases, $x$ is the
+coordinate normal to the interface, $a_2$ is the position of the
+interface and $a_3$ is the width of the
+interface~\cite{Steen-Saethre2014a}.
+This is implemented in {\tt gmx densorder}.
+
+\subsection{Transverse current autocorrelation function}
+In order to establish the transverse current autocorrelation function
+(useful for computing viscosity~\cite{Palmer1994a})
+the following function is fitted:
+\beq
+f(x) ~=~ e^{-\nu}\left({\rm cosh}(\omega\nu)+\frac{{\rm
+ sinh}(\omega\nu)}{\omega}\right)
+\eeq
+with $\nu = x/(2a_0)$ and $\omega = \sqrt{1-a_1}$.
+This is implemented in {\tt gmx tcaf}.
+
+\subsection{Viscosity estimation from pressure autocorrelation
+ function}
+The viscosity is a notoriously difficult property to extract from
+simulations~\cite{Hess2002a,Wensink2003a}. It is {\em in principle}
+possible to determine it by integrating the pressure autocorrelation
+function~\cite{PSmith93c}, however this is often hampered by the noisy
+tail of the ACF. A workaround to this is fitting the ACF to the
+following function~\cite{Guo2002b}:
+\beq
+f(t)/f(0) = (1-C) {\rm cos}(\omega t) e^{-(t/\tau_f)^{\beta_f}} + C
+e^{-(t/\tau_s)^{\beta_s}}
+\eeq
+where $\omega$ is the frequency of rapid pressure oscillations (mainly
+due to bonded forces in molecular simulations), $\tau_f$ and $\beta_f$
+are the time constant and exponent of fast relaxation in a
+stretched-exponential approximation, $\tau_s$ and $\beta_s$ are constants
+for slow relaxation and $C$ is the pre-factor that determines the
+weight between fast and slow relaxation. After a fit, the integral of
+the function $f(t)$ is used to compute the viscosity:
+\beq
+\eta = \frac{V}{k_B T}\int_0^{\infty} f(t) dt
+\eeq
+This equation has been
+applied to computing the bulk and shear viscosity using different
+elements from the pressure tensor~\cite{Fanourgakis2012a}.
+This is implemented in {\tt gmx viscosity}.
+
\section{Mean Square Displacement}
\label{sec:msd}
{\tt gmx msd}\\
keywords = {CHARMM, force field, molecular dynamics, parametrization, DNA, RNA},
year = {2000},
}
+
+@Misc{lmfit,
+ author = {Joachim Wuttke},
+ title = {lmfit},
+ howpublished = {http://apps.jcns.fz-juelich.de/doku/sc/lmfit},
+ year = 2013}
+
+@Article{Steen-Saethre2014a,
+ author = {B. Steen-S{\ae}thre and A. C. Hoffmann and D. van der Spoel},
+ title = {Order Parameters and Algorithmic Approaches for Detection and Demarcation of Interfaces in Hydrate−Fluid and Ice−Fluid Systems},
+ journal = {J. Chem. Theor. Comput.},
+ year = 2014,
+ volume = 10,
+ pages = {5606-5615}}
+
+@Article{ Wensink2003a,
+ author = "E. J. W. Wensink and A. C. Hoffmann and P. J. {van
+ Maaren} and D. {van der Spoel}",
+ title = "Dynamic properties of water/alcohol mixtures studied
+ by computer simulation",
+ journal = BTjcp,
+ year = 2003,
+ volume = 119,
+ pages = "7308-7317"
+}
+
+@Article{Guo2002b,
+ author = {Guang-Jun Guo and Yi-Gang Zhang and Keith Refson and Ya-Juan Zhao},
+ title = {Viscosity and stress autocorrelation function in supercooled water: a molecular dynamics study},
+ journal = {Mol. Phys.},
+ year = 2002,
+ volume = 100,
+ pages = {2617-2627}}
+
+@Article{Fanourgakis2012a,
+ author = {George S. Fanourgakis and J. S. Medina and R. Prosmiti},
+ title = {Determining the Bulk Viscosity of Rigid Water Models},
+ journal = {J. Phys. Chem. A},
+ year = 2012,
+ volume = 116,
+ pages = {2564-2570}}
+
+@Article{Palmer1994a,
+ author = {B. J. Palmer},
+ title = {Transverse-current autocorrelation-function calculations of the shear viscosity for molecular liquids.},
+ journal = {Phys. Rev. E},
+ year = 1994,
+ volume = 49,
+ pages = {359-366}}
+
typedef struct {
const double *t;
const double *y;
+ const double *dy;
double (*f)(double t, const double *par);
} lmcurve_data_struct;
void lmcurve_evaluate( const double *par, int m_dat, const void *data,
double *fvec, gmx_unused int *info )
{
- int i;
+ int i;
+ double fy;
+ lmcurve_data_struct *d = (lmcurve_data_struct*) data;
for (i = 0; i < m_dat; i++)
{
- fvec[i] =
- ((lmcurve_data_struct*)data)->y[i] -
- ((lmcurve_data_struct*)data)->f(
- ((lmcurve_data_struct*)data)->t[i], par );
+ fy = d->f(d->t[i], par );
+ fvec[i] = (d->y[i] - fy)/d->dy[i];
}
}
void lmcurve( int n_par, double *par, int m_dat,
- const double *t, const double *y,
+ const double *t, const double *y, const double *dy,
double (*f)( double t, const double *par ),
const lm_control_struct *control,
lm_status_struct *status )
{
lmcurve_data_struct data;
- data.t = t;
- data.y = y;
- data.f = f;
+ data.t = t;
+ data.y = y;
+ data.dy = dy;
+ data.f = f;
lmmin( n_par, par, m_dat, (const void*) &data,
lmcurve_evaluate, control, status );
__BEGIN_DECLS
void lmcurve( int n_par, double *par, int m_dat,
- const double *t, const double *y,
+ const double *t, const double *y, const double *dy,
double (*f)( double t, const double *par ),
const lm_control_struct *control,
lm_status_struct *status );
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, 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 routine for fitting a data set to a curve
- *
- * \author David van der Spoel <david.vanderspoel@icm.uu.se>
- * \ingroup module_correlationfunctions
- */
-#include "gmxpre.h"
-
-#include "expfit.h"
-
-#include <math.h>
-#include <string.h>
-
-#include "external/lmfit/lmcurve.h"
-
-#include "gromacs/correlationfunctions/integrate.h"
-#include "gromacs/fileio/xvgr.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/real.h"
-#include "gromacs/utility/smalloc.h"
-
-/*! \brief Number of parameters for each fitting function */
-static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3, 6};
-
-const char *s_ffn[effnNR+2] = {
- NULL, "none", "exp", "aexp", "exp_exp", "vac",
- "exp5", "exp7", "exp9", "erffit", "effnPRES", NULL, NULL
-};
-/* We don't allow errest as a choice on the command line */
-
-/*! \brief Long description for each fitting function type */
-static const char *longs_ffn[effnNR] = {
- "no fit",
- "y = exp(-x/a1)",
- "y = a2 exp(-x/a1)",
- "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
- "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
- "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5",
- "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
- "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
- "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
- "y = a2 *2*a1*((exp(-x/a1)-1)*(a1/x)+1)+(1-a2)*2*a3*((exp(-x/a3)-1)*(a3/x)+1)",
- "y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6)"
-};
-
-int effnNparams(int effn)
-{
- if ((0 <= effn) && (effn < effnNR))
- {
- return nfp_ffn[effn];
- }
- else
- {
- return -1;
- }
-}
-
-const char *effnDescription(int effn)
-{
- if ((0 <= effn) && (effn < effnNR))
- {
- return longs_ffn[effn];
- }
- else
- {
- return NULL;
- }
-}
-
-int sffn2effn(const char **sffn)
-{
- int eFitFn, i;
-
- eFitFn = 0;
- for (i = 0; i < effnNR; i++)
- {
- if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0)
- {
- eFitFn = i;
- }
- }
-
- return eFitFn;
-}
-
-/*! \brief Compute exponential function A exp(-x/tau) */
-static double myexp(double x, double A, double tau)
-{
- if ((A == 0) || (tau == 0))
- {
- return 0;
- }
- return A*exp(-x/tau);
-}
-
-/*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
-static double lmc_erffit (double x, const double *a)
-{
- double erfarg;
-
- erfarg = (x-a[2])/(a[3]*a[3]);
- return (a[0]+a[1])/2-(a[0]-a[1])*gmx_erf(erfarg);
-}
-
-/*! \brief Compute y = exp(-x/a0) */
-static double lmc_exp_one_parm(double x, const double *a)
-{
- return exp(-x/a[0]);
-}
-
-/*! \brief Compute y = a1 exp(-x/a0) */
-static double lmc_exp_two_parm(double x, const double *a)
-{
- return a[1]*exp(-x/a[0]);
-}
-
-/*! \brief Compute y = a1 exp(-x/a0) + (1-a1) exp(-x/a2) */
-static double lmc_exp_3_parm(double x, const double *a)
-{
- double e1, e2;
-
- e1 = exp(-x/a[0]);
- e2 = exp(-x/a[2]);
- return a[1]*e1 + (1-a[1])*e2;
-}
-
-/*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 */
-static double lmc_exp_5_parm(double x, const double *a)
-{
- double e1, e2;
-
- e1 = exp(-x/a[1]);
- e2 = exp(-x/a[3]);
- return a[0]*e1 + a[2]*e2 + a[4];
-}
-
-/*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 */
-static double lmc_exp_7_parm(double x, const double *a)
-{
- double e1, e2, e3;
-
- e1 = exp(-x/a[1]);
- e2 = exp(-x/a[3]);
- e3 = exp(-x/a[5]);
- return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6];
-}
-
-/*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 exp(-x/a7) + a8 */
-static double lmc_exp_9_parm(double x, const double *a)
-{
- double e1, e2, e3, e4;
-
- e1 = exp(-x/a[1]);
- e2 = exp(-x/a[3]);
- e3 = exp(-x/a[5]);
- e4 = exp(-x/a[7]);
- return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8];
-}
-
-/*! \brief Compute y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6) */
-static double effnPRES_fun(double x, const double *a)
-{
- double term1, term2, term3;
-
- if (a[0] != 0)
- {
- term3 = a[0] * exp(-pow((x/a[4]), a[5]));
- }
- else
- {
- term3 = 0;
- }
-
- term1 = 1-a[0];
- if (term1 != 0)
- {
- term2 = exp(-pow((x/a[2]), a[3])) * cos(x*a[1]);
- }
- else
- {
- term2 = 0;
- }
-
- return term1*term2 + term3;
-}
-
-/*! \brief Compute vac function */
-static double lmc_vac_2_parm(double x, const double *a)
-{
- /* Fit to function
- *
- * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
- *
- * = exp(-v) (cosh(wv) + 1/w sinh(wv))
- *
- * v = x/(2 a1)
- * w = sqrt(1 - a2)
- *
- * For tranverse current autocorrelation functions:
- * a1 = tau
- * a2 = 4 tau (eta/rho) k^2
- *
- */
-
- double y, v, det, omega, omega2, em, ec, es;
-
- v = x/(2*a[0]);
- det = 1 - a[1];
- em = exp(-v);
- if (det != 0)
- {
- omega2 = fabs(det);
- omega = sqrt(omega2);
- if (det > 0)
- {
- ec = em*0.5*(exp(omega*v)+exp(-omega*v));
- es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
- }
- else
- {
- ec = em*cos(omega*v);
- es = em*sin(omega*v)/omega;
- }
- y = ec + es;
- }
- else
- {
- y = (1+v)*em;
- }
- return y;
-}
-
-/*! \brief Compute error estimate */
-static double lmc_errest_3_parm(double x, const double *a)
-{
- /*
- e1 = (exp(-x/a1) - 1)
- e2 = (exp(-x/a3) - 1)
- v1= 2*a1 * (e1*a1/x + 1)
- v2 = 2*a3 * (e2*a3/x + 1)
- fun = a2*v1 + (1 - a2) * v2
- */
- double e1, e2, v1, v2;
-
- if (a[0] != 0)
- {
- e1 = gmx_expm1(-x/a[0]);
- }
- else
- {
- e1 = 0;
- }
- if (a[2] != 0)
- {
- e2 = gmx_expm1(-x/a[2]);
- }
- else
- {
- e2 = 0;
- }
-
- if (x > 0)
- {
- v1 = 2*a[0]*(e1*a[0]/x + 1);
- v2 = 2*a[2]*(e2*a[2]/x + 1);
- return a[1]*v1 + (1-a[1])*v2;
- }
- else
- {
- return 0;
- }
-}
-
-/*! \brief function type for passing to fitting routine */
-typedef double (*t_lmcurve)(double x, const double *a);
-
-/*! \brief array of fitting functions corresponding to the pre-defined types */
-t_lmcurve lmcurves[effnNR] = {
- lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
- lmc_exp_3_parm, lmc_vac_2_parm,
- lmc_exp_5_parm, lmc_exp_7_parm,
- lmc_exp_9_parm, lmc_erffit, lmc_errest_3_parm, effnPRES_fun
-};
-
-double fit_function(int eFitFn, double *parm, double x)
-{
- // TODO: check that eFitFn is within bounds
- return lmcurves[eFitFn](x, parm);
-}
-
-/*! \brief lmfit_exp supports up to 3 parameter fitting of exponential functions */
-static gmx_bool lmfit_exp(int nfit, double x[], double y[],
- real ftol, int maxiter,
- double parm[], gmx_bool bVerbose,
- int eFitFn)
-{
- double chisq, ochisq;
- gmx_bool bCont;
- int i, j;
- lm_control_struct *control;
- lm_status_struct *status;
-
- if ((eFitFn < 0) || (eFitFn >= effnNR))
- {
- gmx_fatal(FARGS, "fitfn = %d, should be in 0..%d (%s,%d)",
- effnNR-1, eFitFn, __FILE__, __LINE__);
- }
-
- snew(control, 1);
- control->ftol = ftol;
- control->xtol = ftol;
- control->gtol = ftol;
- control->epsilon = 0.1;
- control->stepbound = 100;
- control->patience = maxiter;
- control->scale_diag = 1;
- control->msgfile = stdout;
- control->verbosity = (bVerbose ? 1 : 0);
- control->n_maxpri = nfp_ffn[eFitFn];
- control->m_maxpri = 0;
- snew(status, 1);
- /* Initial params */
- chisq = 1e12;
- j = 0;
- if (bVerbose)
- {
- fprintf(stderr, "%4s %10s Parameters\n",
- "Step", "chi^2");
- }
- do
- {
- ochisq = chisq;
-
- lmcurve(nfp_ffn[eFitFn], parm, nfit, x, y,
- lmcurves[eFitFn], control, status);
- chisq = sqr(status->fnorm);
- if (bVerbose)
- {
- int mmm;
- fprintf(stderr, "%4d %10.5e", j, chisq);
- for (mmm = 0; (mmm < nfp_ffn[eFitFn]); mmm++)
- {
- fprintf(stderr, " %10.5e", parm[mmm]);
- }
- fprintf(stderr, "\n");
- }
- j++;
- bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
- ((ochisq == chisq)));
- }
- while (bCont && (j < maxiter));
- if (bVerbose)
- {
- fprintf(stderr, "\n");
- }
-
- sfree(control);
- sfree(status);
-
- return TRUE;
-}
-
-/*! \brief See description in header file. */
-real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[],
- real begintimefit, real endtimefit, const output_env_t oenv,
- gmx_bool bVerbose, int eFitFn, double fitparms[], int fix)
-{
- FILE *fp;
- char buf[32];
-
- int i, j, nparm, nfitpnts;
- double integral, ttt;
- double *parm, *dparm;
- double *x, *y, *dy;
-#ifdef GMX_DOUBLE
- real ftol = 1e-16;
-#else
- real ftol = 1e-8;
-#endif
- int maxiter = 50;
-
- if (0 != fix)
- {
- fprintf(stderr, "Using fixed parameters in curve fitting is not working anymore\n");
- }
- nparm = nfp_ffn[eFitFn];
- if (debug)
- {
- fprintf(debug, "There are %d points to fit %d vars!\n", ndata, nparm);
- fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
- eFitFn, begintimefit, endtimefit, dt);
- }
-
- snew(x, ndata);
- snew(y, ndata);
- snew(dy, ndata);
-
- j = 0;
- for (i = 0; i < ndata; i++)
- {
- ttt = x0 ? x0[i] : dt*i;
- if (ttt >= begintimefit && ttt <= endtimefit)
- {
- x[j] = ttt;
- y[j] = c1[i];
-
- /* mrqmin does not like sig to be zero */
- if (sig[i] < 1.0e-7)
- {
- dy[j] = 1.0e-7;
- }
- else
- {
- dy[j] = sig[i];
- }
- if (debug)
- {
- fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g\n",
- j, i, x[j], y[j], dy[j]);
- }
- j++;
- }
- }
- nfitpnts = j;
- integral = 0;
- if (nfitpnts < nparm)
- {
- fprintf(stderr, "Not enough data points for fitting!\n");
- }
- else
- {
- snew(parm, nparm);
- snew(dparm, nparm);
- if (fitparms)
- {
- for (i = 0; (i < nparm); i++)
- {
- parm[i] = fitparms[i];
- }
- }
-
- if (!lmfit_exp(nfitpnts, x, y, ftol, maxiter, parm, bVerbose, eFitFn))
- {
- fprintf(stderr, "Fit failed!\n");
- }
- else if (nparm <= 3)
- {
- /* Compute the integral from begintimefit */
- if (nparm == 3)
- {
- integral = (parm[0]*myexp(begintimefit, parm[1], parm[0]) +
- parm[2]*myexp(begintimefit, 1-parm[1], parm[2]));
- }
- else if (nparm == 2)
- {
- integral = parm[0]*myexp(begintimefit, parm[1], parm[0]);
- }
- else if (nparm == 1)
- {
- integral = parm[0]*myexp(begintimefit, 1, parm[0]);
- }
- else
- {
- gmx_fatal(FARGS, "nparm = %d in file %s, line %d",
- nparm, __FILE__, __LINE__);
- }
-
- /* Generate THE output */
- if (bVerbose)
- {
- fprintf(stderr, "FIT: # points used in fit is: %d\n", nfitpnts);
- fprintf(stderr, "FIT: %21s%21s%21s\n",
- "parm0 ", "parm1 (ps) ", "parm2 (ps) ");
- fprintf(stderr, "FIT: ------------------------------------------------------------\n");
- fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
- parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]);
- fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
- begintimefit, integral);
-
- sprintf(buf, "test%d.xvg", nfitpnts);
- fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv);
- fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g\n",
- parm[0], parm[1], parm[2]);
- for (j = 0; (j < nfitpnts); j++)
- {
- ttt = x0 ? x0[j] : dt*j;
- fprintf(fp, "%10.5e %10.5e %10.5e\n",
- ttt, c1[j],
- lmcurves[eFitFn](ttt, parm));
- }
- xvgrclose(fp);
- }
- }
- if (fitparms)
- {
- for (i = 0; (i < nparm); i++)
- {
- fitparms[i] = parm[i];
- }
- }
- sfree(parm);
- sfree(dparm);
- }
-
- sfree(x);
- sfree(y);
- sfree(dy);
-
- return integral;
-}
-
-/*! See description in header file. */
-real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose,
- real tbeginfit, real tendfit, real dt, real c1[], real *fit)
-{
- double fitparm[3];
- double tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate;
- real *sig;
- int i, j, jmax, nf_int;
- gmx_bool bPrint;
-
- bPrint = bVerbose || bDebugMode();
-
- if (bPrint)
- {
- printf("COR:\n");
- }
-
- if (tendfit <= 0)
- {
- tendfit = ncorr*dt;
- }
- nf_int = min(ncorr, (int)(tendfit/dt));
- sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
-
- /* Estimate the correlation time for better fitting */
- ct_estimate = 0.5*c1[0];
- for (i = 1; (i < ncorr) && (c1[i] > 0); i++)
- {
- ct_estimate += c1[i];
- }
- ct_estimate *= dt/c1[0];
-
- if (bPrint)
- {
- printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
- 0.0, dt*nf_int, sum);
- printf("COR: Relaxation times are computed as fit to an exponential:\n");
- printf("COR: %s\n", longs_ffn[fitfn]);
- printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, min(ncorr*dt, tendfit));
- }
-
- tStart = 0;
- if (bPrint)
- {
- printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
- "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
- (nfp_ffn[fitfn] >= 2) ? " a2 ()" : "",
- (nfp_ffn[fitfn] >= 3) ? " a3 (ps)" : "");
- }
-
- snew(sig, ncorr);
-
- if (tbeginfit > 0)
- {
- jmax = 3;
- }
- else
- {
- jmax = 1;
- }
- for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++)
- {
- /* Estimate the correlation time for better fitting */
- c_start = -1;
- ct_estimate = 0;
- for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++)
- {
- if (c_start < 0)
- {
- if (dt*i >= tStart)
- {
- c_start = c1[i];
- ct_estimate = 0.5*c1[i];
- }
- }
- else
- {
- ct_estimate += c1[i];
- }
- }
- if (c_start > 0)
- {
- ct_estimate *= dt/c_start;
- }
- else
- {
- /* The data is strange, so we need to choose somehting */
- ct_estimate = tendfit;
- }
- if (debug)
- {
- fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
- }
-
- if (fitfn == effnEXP3)
- {
- fitparm[0] = 0.002*ncorr*dt;
- fitparm[1] = 0.95;
- fitparm[2] = 0.2*ncorr*dt;
- }
- else
- {
- /* Good initial guess, this increases the probability of convergence */
- fitparm[0] = ct_estimate;
- fitparm[1] = 1.0;
- fitparm[2] = 1.0;
- }
-
- /* Generate more or less appropriate sigma's */
- for (i = 0; i < ncorr; i++)
- {
- sig[i] = sqrt(ct_estimate+dt*i);
- }
-
- nf_int = min(ncorr, (int)((tStart+1e-4)/dt));
- sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
- tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv,
- bDebugMode(), fitfn, fitparm, 0);
- sumtot = sum+tail_corr;
- if (fit && ((jmax == 1) || (j == 1)))
- {
- double mfp[3];
- for (i = 0; (i < 3); i++)
- {
- mfp[i] = fitparm[i];
- }
- for (i = 0; (i < ncorr); i++)
- {
- fit[i] = lmcurves[fitfn](i*dt, mfp);
- }
- }
- if (bPrint)
- {
- printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
- for (i = 0; (i < nfp_ffn[fitfn]); i++)
- {
- printf(" %11.4e", fitparm[i]);
- }
- printf("\n");
- }
- tStart += tbeginfit;
- }
- sfree(sig);
-
- return sumtot;
-}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015, 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 routine for fitting a data set to a curve
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ * \ingroup module_correlationfunctions
+ */
+#include "gmxpre.h"
+
+#include "expfit.h"
+
+#include <math.h>
+#include <string.h>
+
+#include <algorithm>
+
+#include "external/lmfit/lmcurve.h"
+
+#include "gromacs/correlationfunctions/integrate.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/real.h"
+#include "gromacs/utility/smalloc.h"
+
+/*! \brief Number of parameters for each fitting function */
+static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 5, 7, 9, 2, 4, 3, 6 };
+
+/* +2 becuase parse_common_args wants leading and concluding NULL.
+ * We only allow exponential functions as choices on the command line,
+ * hence there are many more NULL field (which have to be at the end of
+ * the array).
+ */
+const char *s_ffn[effnNR+2] = {
+ NULL, "none", "exp", "aexp", "exp_exp",
+ "exp5", "exp7", "exp9",
+ NULL, NULL, NULL, NULL, NULL
+};
+
+/*! \brief Long description for each fitting function type */
+static const char *longs_ffn[effnNR] = {
+ "no fit",
+ "y = exp(-x/|a0|)",
+ "y = a1 exp(-x/|a0|)",
+ "y = a1 exp(-x/|a0|) + (1-a1) exp(-x/(|a2|)), a2 > a0 > 0",
+ "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4, a3 >= a1",
+ "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6, a5 >= a3 >= a1",
+ "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6 exp(-x/|a7|) + a8, a7 >= a5 >= a3 >= a1",
+ "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a0), w = sqrt(1 - a1)",
+ "y = 1/2*(a0+a1) - 1/2*(a0-a1)*erf( (x-a2) /a3^2)",
+ "y = a1 *2*a0*((exp(-x/a0)-1)*(a0/x)+1)+(1-a1)*2*a2*((exp(-x/a2)-1)*(a2/x)+1)",
+ "y = (1-a0)*cos(x*a1)*exp(-(x/a2)^a3) + a0*exp(-(x/a4)^a5)"
+};
+
+int effnNparams(int effn)
+{
+ if ((0 <= effn) && (effn < effnNR))
+ {
+ return nfp_ffn[effn];
+ }
+ else
+ {
+ return -1;
+ }
+}
+
+const char *effnDescription(int effn)
+{
+ if ((0 <= effn) && (effn < effnNR))
+ {
+ return longs_ffn[effn];
+ }
+ else
+ {
+ return NULL;
+ }
+}
+
+int sffn2effn(const char **sffn)
+{
+ int eFitFn, i;
+
+ eFitFn = 0;
+ for (i = 0; i < effnNR; i++)
+ {
+ if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0)
+ {
+ eFitFn = i;
+ }
+ }
+
+ return eFitFn;
+}
+
+/*! \brief Compute exponential function A exp(-x/tau) */
+static double myexp(double x, double A, double tau)
+{
+ if ((A == 0) || (tau == 0))
+ {
+ return 0;
+ }
+ return A*exp(-x/tau);
+}
+
+/*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
+static double lmc_erffit (double x, const double *a)
+{
+ double erfarg;
+ double myerf;
+
+ if (a[3] != 0)
+ {
+ erfarg = (x-a[2])/(a[3]*a[3]);
+ myerf = gmx_erfd(erfarg);
+ }
+ else
+ {
+ /* If a[3] == 0, a[3]^2 = 0 and the erfarg becomes +/- infinity */
+ if (x < a[2])
+ {
+ myerf = -1;
+ }
+ else
+ {
+ myerf = 1;
+ }
+ }
+ return 0.5*((a[0]+a[1]) - (a[0]-a[1])*myerf);
+}
+
+/*! \brief Exponent function that prevents overflow */
+static double safe_exp(double x)
+{
+ double exp_max = 200;
+ double exp_min = -exp_max;
+ if (x <= exp_min)
+ {
+ return exp(exp_min);
+ }
+ else if (x >= exp_max)
+ {
+ return exp(exp_max);
+ }
+ else
+ {
+ return exp(x);
+ }
+}
+
+/*! \brief Exponent minus 1 function that prevents overflow */
+static double safe_expm1(double x)
+{
+ double exp_max = 200;
+ double exp_min = -exp_max;
+ if (x <= exp_min)
+ {
+ return -1;
+ }
+ else if (x >= exp_max)
+ {
+ return exp(exp_max);
+ }
+ else
+ {
+ return gmx_expm1(x);
+ }
+}
+
+/*! \brief Compute y = exp(-x/|a0|) */
+static double lmc_exp_one_parm(double x, const double *a)
+{
+ return safe_exp(-x/fabs(a[0]));
+}
+
+/*! \brief Compute y = a1 exp(-x/|a0|) */
+static double lmc_exp_two_parm(double x, const double *a)
+{
+ return a[1]*safe_exp(-x/fabs(a[0]));
+}
+
+/*! \brief Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|) */
+static double lmc_exp_exp(double x, const double *a)
+{
+ double e1, e2;
+
+ e1 = safe_exp(-x/fabs(a[0]));
+ e2 = safe_exp(-x/(fabs(a[0])+fabs(a[2])));
+ return a[1]*e1 + (1-a[1])*e2;
+}
+
+/*! \brief Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 */
+static double lmc_exp_5_parm(double x, const double *a)
+{
+ double e1, e2;
+
+ e1 = safe_exp(-x/fabs(a[1]));
+ e2 = safe_exp(-x/(fabs(a[1])+fabs(a[3])));
+ return a[0]*e1 + a[2]*e2 + a[4];
+}
+
+/*! \brief Compute 7 parameter exponential function value.
+ *
+ * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
+ * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6
+ */
+static double lmc_exp_7_parm(double x, const double *a)
+{
+ double e1, e2, e3;
+ double fa1, fa3, fa5;
+
+ fa1 = fabs(a[1]);
+ fa3 = fa1 + fabs(a[3]);
+ fa5 = fa3 + fabs(a[5]);
+ e1 = safe_exp(-x/fa1);
+ e2 = safe_exp(-x/fa3);
+ e3 = safe_exp(-x/fa5);
+ return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6];
+}
+
+/*! \brief Compute 9 parameter exponential function value.
+ *
+ * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
+ * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 exp(-x/(|a1|+|a3|+|a5|+|a7|)) + a8
+ */
+static double lmc_exp_9_parm(double x, const double *a)
+{
+ double e1, e2, e3, e4;
+ double fa1, fa3, fa5, fa7;
+
+ fa1 = fabs(a[1]);
+ fa3 = fa1 + fabs(a[3]);
+ fa5 = fa3 + fabs(a[5]);
+ fa7 = fa5 + fabs(a[7]);
+
+ e1 = safe_exp(-x/fa1);
+ e2 = safe_exp(-x/fa3);
+ e3 = safe_exp(-x/fa5);
+ e4 = safe_exp(-x/fa7);
+ return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8];
+}
+
+/*! \brief Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|) */
+static double lmc_pres_6_parm(double x, const double *a)
+{
+ double term1, term2, term3;
+ double pow_max = 10;
+
+ term3 = 0;
+ if ((a[4] != 0) && (a[0] != 0))
+ {
+ double power = std::min(fabs(a[5]), pow_max);
+ term3 = a[0] * safe_exp(-pow((x/fabs(a[4])), power));
+ }
+
+ term1 = 1-a[0];
+ term2 = 0;
+ if ((term1 != 0) && (a[2] != 0))
+ {
+ double power = std::min(fabs(a[3]), pow_max);
+ term2 = safe_exp(-pow((x/fabs(a[2])), power)) * cos(x*fabs(a[1]));
+ }
+
+ return term1*term2 + term3;
+}
+
+/*! \brief Compute vac function */
+static double lmc_vac_2_parm(double x, const double *a)
+{
+ /* Fit to function
+ *
+ * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
+ *
+ * = exp(-v) (cosh(wv) + 1/w sinh(wv))
+ *
+ * v = x/(2 a0)
+ * w = sqrt(1 - a1)
+ *
+ * For tranverse current autocorrelation functions:
+ * a0 = tau
+ * a1 = 4 tau (eta/rho) k^2
+ *
+ */
+
+ double y, v, det, omega, wv, em, ec, es;
+ double wv_max = 100;
+
+ v = x/(2*fabs(a[0]));
+ det = 1 - a[1];
+ em = safe_exp(-v);
+ if (det != 0)
+ {
+ omega = sqrt(fabs(det));
+ wv = std::min(omega*v, wv_max);
+
+ if (det > 0)
+ {
+ ec = em*0.5*(safe_exp(wv)+safe_exp(-wv));
+ es = em*0.5*(safe_exp(wv)-safe_exp(-wv))/omega;
+ }
+ else
+ {
+ ec = em*cos(wv);
+ es = em*sin(wv)/omega;
+ }
+ y = ec + es;
+ }
+ else
+ {
+ y = (1+v)*em;
+ }
+ return y;
+}
+
+/*! \brief Compute error estimate */
+static double lmc_errest_3_parm(double x, const double *a)
+{
+ double e1, e2, v1, v2;
+ double fa0 = fabs(a[0]);
+ double fa1;
+ double fa2 = fa0+fabs(a[2]);
+
+ if (a[0] != 0)
+ {
+ e1 = safe_expm1(-x/fa0);
+ }
+ else
+ {
+ e1 = 0;
+ }
+ if (a[2] != 0)
+ {
+ e2 = safe_expm1(-x/fa2);
+ }
+ else
+ {
+ e2 = 0;
+ }
+
+ if (x > 0)
+ {
+ v1 = 2*fa0*(e1*fa0/x + 1);
+ v2 = 2*fa2*(e2*fa2/x + 1);
+ /* We need 0 <= a1 <= 1 */
+ fa1 = std::min(1.0, std::max(0.0, a[1]));
+
+ return fa1*v1 + (1-fa1)*v2;
+ }
+ else
+ {
+ return 0;
+ }
+}
+
+/*! \brief function type for passing to fitting routine */
+typedef double (*t_lmcurve)(double x, const double *a);
+
+/*! \brief array of fitting functions corresponding to the pre-defined types */
+t_lmcurve lmcurves[effnNR+1] = {
+ lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
+ lmc_exp_exp, lmc_exp_5_parm, lmc_exp_7_parm,
+ lmc_exp_9_parm,
+ lmc_vac_2_parm, lmc_erffit, lmc_errest_3_parm, lmc_pres_6_parm
+};
+
+double fit_function(const int eFitFn, const double parm[], const double x)
+{
+ if ((eFitFn < 0) || (eFitFn >= effnNR))
+ {
+ fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n",
+ eFitFn, effnNR-1);
+ return 0.0;
+ }
+ return lmcurves[eFitFn](x, parm);
+}
+
+/*! \brief lmfit_exp supports fitting of different functions
+ *
+ * This routine calls the Levenberg-Marquardt non-linear fitting
+ * routine for fitting a data set with errors to a target function.
+ * Fitting routines included in gromacs in src/external/lmfit.
+ */
+static gmx_bool lmfit_exp(int nfit,
+ const double x[],
+ const double y[],
+ const double dy[],
+ double parm[],
+ gmx_bool bVerbose,
+ int eFitFn,
+ int nfix)
+{
+ double chisq, ochisq;
+ gmx_bool bCont;
+ int j;
+ int maxiter = 100;
+ lm_control_struct control;
+ lm_status_struct *status;
+ int nparam = effnNparams(eFitFn);
+ int p2;
+ gmx_bool bSkipLast;
+
+ if ((eFitFn < 0) || (eFitFn >= effnNR))
+ {
+ fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n",
+ eFitFn, effnNR-1);
+ return FALSE;
+ }
+ /* Using default control structure for double precision fitting that
+ * comes with the lmfit package (i.e. from the include file).
+ */
+ control = lm_control_double;
+ control.verbosity = (bVerbose ? 1 : 0);
+ control.n_maxpri = 0;
+ control.m_maxpri = 0;
+
+ snew(status, 1);
+ /* Initial params */
+ chisq = 1e12;
+ j = 0;
+ if (bVerbose)
+ {
+ printf("%4s %10s Parameters\n", "Step", "chi^2");
+ }
+ /* Check whether we have to skip some params */
+ if (nfix > 0)
+ {
+ do
+ {
+ p2 = 1 << (nparam-1);
+ bSkipLast = ((p2 & nfix) == p2);
+ if (bSkipLast)
+ {
+ nparam--;
+ nfix -= p2;
+ }
+ }
+ while ((nparam > 0) && (bSkipLast));
+ if (bVerbose)
+ {
+ printf("Using %d out of %d parameters\n", nparam, effnNparams(eFitFn));
+ }
+ }
+ do
+ {
+ ochisq = chisq;
+ lmcurve(nparam, parm, nfit, x, y, dy,
+ lmcurves[eFitFn], &control, status);
+ chisq = sqr(status->fnorm);
+ if (bVerbose)
+ {
+ printf("status: fnorm = %g, nfev = %d, userbreak = %d\noutcome = %s\n",
+ status->fnorm, status->nfev, status->userbreak,
+ lm_infmsg[status->outcome]);
+ }
+ if (bVerbose)
+ {
+ int mmm;
+ printf("%4d %8g", j, chisq);
+ for (mmm = 0; (mmm < effnNparams(eFitFn)); mmm++)
+ {
+ printf(" %8g", parm[mmm]);
+ }
+ printf("\n");
+ }
+ j++;
+ bCont = (fabs(ochisq - chisq) > fabs(control.ftol*chisq));
+ }
+ while (bCont && (j < maxiter));
+
+ sfree(status);
+
+ return TRUE;
+}
+
+/*! \brief Ensure the fitting parameters are well-behaved.
+ *
+ * In order to make sure that fitting behaves according to the
+ * constraint that time constants are positive and increasing
+ * we enforce this here by setting all time constants to their
+ * absolute value and by adding e.g. |a_0| to |a_2|. This is
+ * done in conjunction with the fitting functions themselves.
+ * When there are multiple time constants we make sure that
+ * the successive time constants are at least double the
+ * previous ones and during fitting we enforce the they remain larger.
+ * It may very well help the convergence of the fitting routine.
+ */
+static void initiate_fit_params(int eFitFn,
+ double params[])
+{
+ int i, nparm;
+
+ nparm = effnNparams(eFitFn);
+
+ switch (eFitFn)
+ {
+ case effnVAC:
+ GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
+ break;
+ case effnEXP1:
+ case effnEXP2:
+ case effnEXPEXP:
+ GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
+ if (nparm > 2)
+ {
+ GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
+ /* In order to maintain params[2] >= params[0] in the final
+ * result, we fit the difference between the two, that is
+ * params[2]-params[0] and in the add add in params[0]
+ * again.
+ */
+ params[2] = std::max(fabs(params[2])-params[0], params[0]);
+ }
+ break;
+ case effnEXP5:
+ case effnEXP7:
+ case effnEXP9:
+ GMX_ASSERT(params[1] >= 0, "parameters should be >= 0");
+ params[1] = fabs(params[1]);
+ if (nparm > 3)
+ {
+ GMX_ASSERT(params[3] >= 0, "parameters should be >= 0");
+ /* See comment under effnEXPEXP */
+ params[3] = std::max(fabs(params[3])-params[1], params[1]);
+ if (nparm > 5)
+ {
+ GMX_ASSERT(params[5] >= 0, "parameters should be >= 0");
+ /* See comment under effnEXPEXP */
+ params[5] = std::max(fabs(params[5])-params[3], params[3]);
+ if (nparm > 7)
+ {
+ GMX_ASSERT(params[7] >= 0, "parameters should be >= 0");
+ /* See comment under effnEXPEXP */
+ params[7] = std::max(fabs(params[7])-params[5], params[5]);
+ }
+ }
+ }
+ break;
+ case effnERREST:
+ GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
+ GMX_ASSERT(params[1] >= 0 && params[1] <= 1, "parameter 1 should in 0 .. 1");
+ GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
+ /* See comment under effnEXPEXP */
+ params[2] = fabs(params[2])-params[0];
+ break;
+ case effnPRES:
+ for (i = 1; (i < nparm); i++)
+ {
+ GMX_ASSERT(params[i] >= 0, "parameters should be >= 0");
+ }
+ break;
+ default:
+ break;
+ }
+}
+
+/*! \brief Process the fitting parameters to get output parameters.
+ *
+ * See comment at the previous function.
+ */
+static void extract_fit_params(int eFitFn,
+ double params[])
+{
+ int i, nparm;
+
+ nparm = effnNparams(eFitFn);
+
+ switch (eFitFn)
+ {
+ case effnVAC:
+ params[0] = fabs(params[0]);
+ break;
+ case effnEXP1:
+ case effnEXP2:
+ case effnEXPEXP:
+ params[0] = fabs(params[0]);
+ if (nparm > 2)
+ {
+ /* Back conversion of parameters from the fitted difference
+ * to the absolute value.
+ */
+ params[2] = fabs(params[2])+params[0];
+ }
+ break;
+ case effnEXP5:
+ case effnEXP7:
+ case effnEXP9:
+ params[1] = fabs(params[1]);
+ if (nparm > 3)
+ {
+ /* See comment under effnEXPEXP */
+ params[3] = fabs(params[3])+params[1];
+ if (nparm > 5)
+ {
+ /* See comment under effnEXPEXP */
+ params[5] = fabs(params[5])+params[3];
+ if (nparm > 7)
+ {
+ /* See comment under effnEXPEXP */
+ params[7] = fabs(params[7])+params[5];
+ }
+ }
+ }
+ break;
+ case effnERREST:
+ params[0] = fabs(params[0]);
+ if (params[1] < 0)
+ {
+ params[1] = 0;
+ }
+ else if (params[1] > 1)
+ {
+ params[1] = 1;
+ }
+ /* See comment under effnEXPEXP */
+ params[2] = params[0]+fabs(params[2]);
+ break;
+ case effnPRES:
+ for (i = 1; (i < nparm); i++)
+ {
+ params[i] = fabs(params[i]);
+ }
+ break;
+ default:
+ break;
+ }
+}
+
+/*! \brief Print chi-squared value and the parameters */
+static void print_chi2_params(FILE *fp,
+ const int eFitFn,
+ const double fitparms[],
+ const char *label,
+ const int nfitpnts,
+ const double x[],
+ const double y[])
+{
+ int i;
+ double chi2 = 0;
+
+ for (i = 0; (i < nfitpnts); i++)
+ {
+ double yfit = lmcurves[eFitFn](x[i], fitparms);
+ chi2 += sqr(y[i] - yfit);
+ }
+ fprintf(fp, "There are %d data points, %d parameters, %s chi2 = %g\nparams:",
+ nfitpnts, effnNparams(eFitFn), label, chi2);
+ for (i = 0; (i < effnNparams(eFitFn)); i++)
+ {
+ fprintf(fp, " %10g", fitparms[i]);
+ }
+ fprintf(fp, "\n");
+}
+
+/*! \brief See description in header file. */
+real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[],
+ real begintimefit, real endtimefit, const output_env_t oenv,
+ gmx_bool bVerbose, int eFitFn, double fitparms[], int fix,
+ const char *fn_fitted)
+{
+ FILE *fp;
+ int i, j, nfitpnts;
+ double integral, ttt;
+ double *x, *y, *dy;
+
+ if (0 != fix)
+ {
+ fprintf(stderr, "Using fixed parameters in curve fitting is temporarily not working.\n");
+ }
+ if (debug)
+ {
+ fprintf(debug, "There are %d points to fit %d vars!\n", ndata, effnNparams(eFitFn));
+ fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
+ eFitFn, begintimefit, endtimefit, dt);
+ }
+
+ snew(x, ndata);
+ snew(y, ndata);
+ snew(dy, ndata);
+ j = 0;
+ for (i = 0; i < ndata; i++)
+ {
+ ttt = x0 ? x0[i] : dt*i;
+ if (ttt >= begintimefit && ttt <= endtimefit)
+ {
+ x[j] = ttt;
+ y[j] = c1[i];
+ if (NULL == sig)
+ {
+ // No weighting if all values are divided by 1.
+ dy[j] = 1;
+ }
+ else
+ {
+ dy[j] = std::max(1.0e-7, (double)sig[i]);
+ }
+ if (debug)
+ {
+ fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy=%g, ttt=%g\n",
+ j, i, x[j], y[j], dy[j], ttt);
+ }
+ j++;
+ }
+ }
+ nfitpnts = j;
+ integral = 0;
+ if (nfitpnts < effnNparams(eFitFn))
+ {
+ fprintf(stderr, "Not enough (%d) data points for fitting, dt = %g!\n",
+ nfitpnts, dt);
+ }
+ else
+ {
+ gmx_bool bSuccess;
+
+ if (bVerbose)
+ {
+ print_chi2_params(stdout, eFitFn, fitparms, "initial", nfitpnts, x, y);
+ }
+ initiate_fit_params(eFitFn, fitparms);
+
+ bSuccess = lmfit_exp(nfitpnts, x, y, dy, fitparms, bVerbose, eFitFn, fix);
+ extract_fit_params(eFitFn, fitparms);
+
+ if (!bSuccess)
+ {
+ fprintf(stderr, "Fit failed!\n");
+ }
+ else
+ {
+ if (bVerbose)
+ {
+ print_chi2_params(stdout, eFitFn, fitparms, "final", nfitpnts, x, y);
+ }
+ switch (eFitFn)
+ {
+ case effnEXP1:
+ integral = fitparms[0]*myexp(begintimefit, 1, fitparms[0]);
+ break;
+ case effnEXP2:
+ integral = fitparms[0]*myexp(begintimefit, fitparms[1], fitparms[0]);
+ break;
+ case effnEXPEXP:
+ integral = (fitparms[0]*myexp(begintimefit, fitparms[1], fitparms[0]) +
+ fitparms[2]*myexp(begintimefit, 1-fitparms[1], fitparms[2]));
+ break;
+ case effnEXP5:
+ case effnEXP7:
+ case effnEXP9:
+ integral = 0;
+ for (i = 0; (i < (effnNparams(eFitFn)-1)/2); i++)
+ {
+ integral += fitparms[2*i]*myexp(begintimefit, fitparms[2*i+1], fitparms[2*i]);
+ }
+ break;
+ default:
+ /* Do numerical integration */
+ integral = 0;
+ for (i = 0; (i < nfitpnts-1); i++)
+ {
+ double y0 = lmcurves[eFitFn](x[i], fitparms);
+ double y1 = lmcurves[eFitFn](x[i+1], fitparms);
+ integral += (x[i+1]-x[i])*(y1+y0)*0.5;
+ }
+ break;
+ }
+
+ if (bVerbose)
+ {
+ printf("FIT: Integral of fitted function: %g\n", integral);
+ if ((effnEXP5 == eFitFn) || (effnEXP7 == eFitFn) || (effnEXP9 == eFitFn))
+ {
+ printf("FIT: Note that the constant term is not taken into account when computing integral.\n");
+ }
+ }
+ /* Generate debug output */
+ if (NULL != fn_fitted)
+ {
+ fp = xvgropen(fn_fitted, "Data + Fit", "Time (ps)",
+ "Data (t)", oenv);
+ for (i = 0; (i < effnNparams(eFitFn)); i++)
+ {
+ fprintf(fp, "# fitparms[%d] = %g\n", i, fitparms[i]);
+ }
+ for (j = 0; (j < nfitpnts); j++)
+ {
+ real ttt = x0 ? x0[i] : dt*j;
+ fprintf(fp, "%10.5e %10.5e %10.5e\n",
+ x[j], y[j], lmcurves[eFitFn](ttt, fitparms));
+ }
+ xvgrclose(fp);
+ }
+ }
+ }
+
+ sfree(x);
+ sfree(y);
+ sfree(dy);
+
+ return integral;
+}
+
+/*! See description in header file. */
+real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose,
+ real tbeginfit, real tendfit, real dt, real c1[], real *fit)
+{
+ double fitparm[3];
+ double tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate;
+ real *sig;
+ int i, j, jmax, nf_int;
+ gmx_bool bPrint;
+
+ bPrint = bVerbose || bDebugMode();
+
+ if (bPrint)
+ {
+ printf("COR:\n");
+ }
+
+ if (tendfit <= 0)
+ {
+ tendfit = ncorr*dt;
+ }
+ nf_int = std::min(ncorr, (int)(tendfit/dt));
+ sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
+
+ /* Estimate the correlation time for better fitting */
+ ct_estimate = 0.5*c1[0];
+ for (i = 1; (i < ncorr) && (c1[i] > 0); i++)
+ {
+ ct_estimate += c1[i];
+ }
+ ct_estimate *= dt/c1[0];
+
+ if (bPrint)
+ {
+ printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
+ 0.0, dt*nf_int, sum);
+ printf("COR: Relaxation times are computed as fit to an exponential:\n");
+ printf("COR: %s\n", effnDescription(fitfn));
+ printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, std::min(ncorr*dt, tendfit));
+ }
+
+ tStart = 0;
+ if (bPrint)
+ {
+ printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
+ "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
+ (effnNparams(fitfn) >= 2) ? " a2 ()" : "",
+ (effnNparams(fitfn) >= 3) ? " a3 (ps)" : "");
+ }
+
+ snew(sig, ncorr);
+
+ if (tbeginfit > 0)
+ {
+ jmax = 3;
+ }
+ else
+ {
+ jmax = 1;
+ }
+ for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++)
+ {
+ /* Estimate the correlation time for better fitting */
+ c_start = -1;
+ ct_estimate = 0;
+ for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++)
+ {
+ if (c_start < 0)
+ {
+ if (dt*i >= tStart)
+ {
+ c_start = c1[i];
+ ct_estimate = 0.5*c1[i];
+ }
+ }
+ else
+ {
+ ct_estimate += c1[i];
+ }
+ }
+ if (c_start > 0)
+ {
+ ct_estimate *= dt/c_start;
+ }
+ else
+ {
+ /* The data is strange, so we need to choose somehting */
+ ct_estimate = tendfit;
+ }
+ if (debug)
+ {
+ fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
+ }
+
+ if (fitfn == effnEXPEXP)
+ {
+ fitparm[0] = 0.002*ncorr*dt;
+ fitparm[1] = 0.95;
+ fitparm[2] = 0.2*ncorr*dt;
+ }
+ else
+ {
+ /* Good initial guess, this increases the probability of convergence */
+ fitparm[0] = ct_estimate;
+ fitparm[1] = 1.0;
+ fitparm[2] = 1.0;
+ }
+
+ /* Generate more or less appropriate sigma's */
+ for (i = 0; i < ncorr; i++)
+ {
+ sig[i] = sqrt(ct_estimate+dt*i);
+ }
+
+ nf_int = std::min(ncorr, (int)((tStart+1e-4)/dt));
+ sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
+ tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv,
+ bDebugMode(), fitfn, fitparm, 0, NULL);
+ sumtot = sum+tail_corr;
+ if (fit && ((jmax == 1) || (j == 1)))
+ {
+ double mfp[3];
+ for (i = 0; (i < 3); i++)
+ {
+ mfp[i] = fitparm[i];
+ }
+ for (i = 0; (i < ncorr); i++)
+ {
+ fit[i] = lmcurves[fitfn](i*dt, mfp);
+ }
+ }
+ if (bPrint)
+ {
+ printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
+ for (i = 0; (i < effnNparams(fitfn)); i++)
+ {
+ printf(" %11.4e", fitparm[i]);
+ }
+ printf("\n");
+ }
+ tStart += tbeginfit;
+ }
+ sfree(sig);
+
+ return sumtot;
+}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015, 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.
* Enum to select fitting functions
*/
enum {
- effnNONE, effnEXP1, effnEXP2, effnEXP3, effnVAC,
- effnEXP5, effnEXP7, effnEXP9, effnERF, effnERREST, effnPRES, effnNR
+ effnNONE, effnEXP1, effnEXP2, effnEXPEXP,
+ effnEXP5, effnEXP7, effnEXP9,
+ effnVAC, effnERF, effnERREST, effnPRES, effnNR
};
/*! \brief
* \param[in] x The value of x
* \return the value of the fit
*/
-double fit_function(int eFitFn, double *parm, double x);
+double fit_function(const int eFitFn, const double parm[], const double x);
/*! \brief
* Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential
* If x == NULL, the timestep dt will be used to create a time axis.
* \param[in] ndata Number of data points
* \param[in] c1 The data points
- * \param[in] sig The standard deviation in the points
+ * \param[in] sig The standard deviation in the points (can be NULL)
* \param[in] dt The time step
* \param[in] x The X-axis (may be NULL, see above)
* \param[in] begintimefit Starting time for fitting
* \param[in] oenv Output formatting information
* \param[in] bVerbose Should the routine write to console?
* \param[in] eFitFn Fitting function (0 .. effnNR)
- * \param[out] fitparms[]
+ * \param[inout] fitparms[] Fitting parameters, see printed manual for a
+ * detailed description. Note that in all implemented cases the parameters
+ * corresponding to time constants will be generated with increasing values.
+ * Such input parameters should therefore be provided in increasing order.
+ * If this is not the case or if subsequent time parameters differ by less than
+ * a factor of 2, they will be modified to ensure tau_i+1 >= 2 tau_i.
* \param[in] fix Constrains fit parameter i at it's starting value, when the i'th bit
- * of fix is set.
+ * of fix is set. This works only when the N last parameters are fixed
+ * but not when a parameter somewhere in the middle needs to be fixed.
+ * \param[in] fn_fitted If not NULL file to print the data and fitted curve to
* \return integral.
*/
real do_lmfit(int ndata, real c1[], real sig[], real dt, real *x,
real begintimefit, real endtimefit, const output_env_t oenv,
- gmx_bool bVerbose, int eFitFn, double fitparms[], int fix);
+ gmx_bool bVerbose, int eFitFn, double fitparms[], int fix,
+ const char *fn_fitted);
/*! \brief
* Fit an autocorrelation function to a pre-defined functional form
* \param[in] tendfit Ending time for fitting
* \param[in] dt The time step
* \param[in] c1 The data points
- * \param[out] fit The fitting parameters
+ * \param[inout] fit The fitting parameters
* \return the integral over the autocorrelation function?
*/
real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose,
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015, 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.
#include <gtest/gtest.h>
#include "gromacs/fileio/xvgr.h"
+#include "gromacs/legacyheaders/oenv.h"
#include "gromacs/utility/smalloc.h"
#include "testutils/refdata.h"
#include "testutils/testasserts.h"
#include "testutils/testfilemanager.h"
-//! Number of data files for testing.
-#define expTestNrTypes 3
-
namespace gmx
{
namespace
{
+class ExpfitData
+{
+ public:
+ int nrLines_;
+ std::vector<real> x_, y_;
+ real startTime_, endTime_, dt_;
+};
+
class ExpfitTest : public ::testing::Test
{
protected:
- static int nrLines_;
- static std::vector<real> values_[expTestNrTypes];
- static int nrColumns_;
- static std::vector<real> standardDev_;
- static real startTime_;
- static real endTime_;
- static real timeDeriv_;
- test::TestReferenceData refData_;
- test::TestReferenceChecker checker_;
+ static std::vector<ExpfitData> data_;
+ test::TestReferenceData refData_;
+ test::TestReferenceChecker checker_;
ExpfitTest( )
: checker_(refData_.rootChecker())
{
// Static initiation, only run once every test.
static void SetUpTestCase()
{
- double ** tempValues;
- std::string fileName[expTestNrTypes];
- fileName[0] = test::TestFileManager::getInputFilePath("testINVEXP.xvg");
- fileName[1] = test::TestFileManager::getInputFilePath("testPRES.xvg");
- fileName[2] = test::TestFileManager::getInputFilePath("testEXP.xvg");
- for (int i = 0; i < expTestNrTypes; i++)
+ double ** tempValues = NULL;
+ std::vector<std::string> fileName;
+ fileName.push_back(test::TestFileManager::getInputFilePath("testINVEXP.xvg"));
+ fileName.push_back(test::TestFileManager::getInputFilePath("testPRES.xvg"));
+ fileName.push_back(test::TestFileManager::getInputFilePath("testINVEXP79.xvg"));
+ fileName.push_back(test::TestFileManager::getInputFilePath("testERF.xvg"));
+ fileName.push_back(test::TestFileManager::getInputFilePath("testERREST.xvg"));
+ for (std::vector<std::string>::iterator i = fileName.begin(); i < fileName.end(); ++i)
{
- const char * name = fileName[i].c_str();
- // TODO: this assumes all files have the same length.
- nrLines_ = read_xvg(name, &tempValues, &nrColumns_);
-
- // Generating standard deviation
- if (i == 0)
+ const char * name = i->c_str();
+ int nrColumns;
+ ExpfitData ed;
+ ed.nrLines_ = read_xvg(name, &tempValues, &nrColumns);
+ ed.dt_ = tempValues[0][1] - tempValues[0][0];
+ ed.startTime_ = tempValues[0][0];
+ ed.endTime_ = tempValues[0][ed.nrLines_-1];
+ for (int j = 0; j < ed.nrLines_; j++)
{
- double fac = 1.0/nrLines_;
- for (int j = 0; j < nrLines_; j++)
- {
- standardDev_.push_back(fac);
- }
- timeDeriv_ = tempValues[0][1] - tempValues[0][0];
- startTime_ = tempValues[0][0];
- endTime_ = tempValues[0][nrLines_-1];
- }
-
- for (int j = 0; j < nrLines_; j++)
- {
- values_[i].push_back((real)tempValues[1][j]);
+ ed.x_.push_back((real)tempValues[0][j]);
+ ed.y_.push_back((real)tempValues[1][j]);
}
+ data_.push_back(ed);
// Free memory that was allocated in read_xvg
- for (int i = 0; i < nrColumns_; i++)
+ for (int j = 0; j < nrColumns; j++)
{
- sfree(tempValues[i]);
- tempValues[i] = NULL;
+ sfree(tempValues[j]);
+ tempValues[j] = NULL;
}
sfree(tempValues);
tempValues = NULL;
{
}
- void test(int type, double result[], double tolerance, int testType)
+ void test(int type, double result[], double tolerance,
+ unsigned int testType)
{
- int nfitparm = effnNparams(type);
-
- do_lmfit(nrLines_, &values_[testType][0], &standardDev_[0], timeDeriv_,
- NULL, startTime_, endTime_, NULL, false, type, result, 0);
+ int nfitparm = effnNparams(type);
+ output_env_t oenv;
+ if (testType >= data_.size())
+ {
+ GMX_THROW(InvalidInputError("testType out of range"));
+ }
+ output_env_init_default(&oenv);
+ do_lmfit(data_[testType].nrLines_,
+ &(data_[testType].y_[0]),
+ NULL,
+ data_[testType].dt_,
+ &(data_[testType].x_[0]),
+ data_[testType].startTime_,
+ data_[testType].endTime_,
+ oenv, false, type, result, 0, NULL);
+ output_env_done(oenv);
checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, tolerance));
checker_.checkSequenceArray(nfitparm, result, "result");
}
//static var
-int ExpfitTest::nrLines_;
-//cppcheck-suppress arrayIndexOutOfBounds fixed in 1.68-dev
-std::vector<real> ExpfitTest::values_[expTestNrTypes];
-int ExpfitTest::nrColumns_;
-std::vector<real> ExpfitTest::standardDev_;
-real ExpfitTest::startTime_;
-real ExpfitTest::endTime_;
-real ExpfitTest::timeDeriv_;
+std::vector<ExpfitData> ExpfitTest::data_;
TEST_F (ExpfitTest, EffnEXP1) {
double param[] = {25};
- test(effnEXP1, param, 1e-6, 0);
+ test(effnEXP1, param, 1e-5, 0);
}
TEST_F (ExpfitTest, EffnEXP2) {
double param[] = {35, 0.5};
- test(effnEXP2, param, 1e-6, 0);
+ test(effnEXP2, param, 3e-5, 0);
}
-TEST_F (ExpfitTest, EffnEXP3) {
- double param[] = {45, 0.5, 5};
- test(effnEXP3, param, 1e-4, 0);
+TEST_F (ExpfitTest, EffnEXPEXP) {
+ double param[] = {5, 0.5, 45};
+ test(effnEXPEXP, param, 1e-2, 0);
}
TEST_F (ExpfitTest, EffnEXP5) {
double param[] = {0.5, 5, 0.5, 50, 0.002};
- test(effnEXP5, param, 1e-4, 0);
+ test(effnEXP5, param, 1e-2, 2);
}
TEST_F (ExpfitTest, EffnEXP7) {
- double param[] = {0.5, 5, -0.02, 0.5, 0.5, 50, -0.002};
- test(effnEXP7, param, 1e-4, 0);
+ double param[] = {0.1, 2, 0.5, 30, 0.3, 50, -0.002};
+ test(effnEXP7, param, 1e-2, 2);
}
TEST_F (ExpfitTest, EffnEXP9) {
- double param[] = {2, 1200, -1, 300, 0.7, 70, 0.5, 6, -0.5};
- test(effnEXP9, param, 4e-2, 0);
+ double param[] = {0.4, 5, 0.2, 30, 0.1, 70, 0.2, 200, -0.05};
+ test(effnEXP9, param, 4e-2, 2);
}
TEST_F (ExpfitTest, EffnERF) {
- double param[] = {0.5, 0.5, 0.5, 1};
- test(effnERF, param, 1e-2, 0);
+ double param[] = {80, 120, 180, 5};
+ test(effnERF, param, 1e-1, 3);
}
TEST_F (ExpfitTest, EffnERREST) {
- double param[] = {0.5, 0.7, 0.3};
- test(effnERREST, param, 1e-4, 2);
+ double param[] = {1, 0.9, 100};
+ test(effnERREST, param, 5e-3, 4);
}
TEST_F (ExpfitTest, EffnVAC) {
- double param[] = {0.5, 0.05};
- test(effnVAC, param, 1e-4, 0);
+ double param[] = {30, 0.0};
+ test(effnVAC, param, 0.05, 0);
}
-TEST_F (ExpfitTest, DISABLED_EffnPRES) {
- //TODO: This test is prodocues NaNs and INFs. Fix and then reactivate.
- double param[] = {0, 10, 4, 1, 0.5, 1};
+TEST_F (ExpfitTest, EffnPRES) {
+ double param[] = {0.6, 10, 7, 1, 0.25, 2};
test(effnPRES, param, 1e-4, 1);
}
<Real>0.98145622</Real>
<Real>0.97775453</Real>
<Real>0.97388524</Real>
- <Real>0.96959651</Real>
+ <Real>0.96959645</Real>
<Real>0.96521682</Real>
<Real>0.96071005</Real>
<Real>0.95656222</Real>
<Real>0.92004192</Real>
<Real>0.92134565</Real>
<Real>0.92249507</Real>
- <Real>0.92415452</Real>
+ <Real>0.92415446</Real>
<Real>0.92598462</Real>
<Real>0.92835307</Real>
<Real>0.93056071</Real>
<Real>0.9445464</Real>
<Real>0.94751155</Real>
<Real>0.95099658</Real>
- <Real>0.95437491</Real>
- <Real>0.95766789</Real>
+ <Real>0.95437485</Real>
+ <Real>0.95766783</Real>
<Real>0.96085715</Real>
<Real>0.96405703</Real>
<Real>0.96722084</Real>
<Real>0.98345906</Real>
<Real>0.98124063</Real>
<Real>0.97880316</Real>
- <Real>0.97668463</Real>
+ <Real>0.97668457</Real>
<Real>0.97459376</Real>
<Real>0.97177458</Real>
<Real>0.9690451</Real>
<Real>0.96629435</Real>
- <Real>0.96378684</Real>
+ <Real>0.96378678</Real>
<Real>0.96081358</Real>
<Real>0.95825708</Real>
<Real>0.95578271</Real>
<Real>0.95307195</Real>
- <Real>0.95074266</Real>
+ <Real>0.9507426</Real>
<Real>0.94830912</Real>
- <Real>0.94601887</Real>
+ <Real>0.94601882</Real>
<Real>0.94421709</Real>
<Real>0.94242907</Real>
<Real>0.94094354</Real>
<Real>0.94294119</Real>
<Real>0.94477904</Real>
<Real>0.94679701</Real>
- <Real>0.94860345</Real>
+ <Real>0.94860339</Real>
<Real>0.95079654</Real>
- <Real>0.95264018</Real>
+ <Real>0.95264012</Real>
<Real>0.95462155</Real>
<Real>0.95684212</Real>
<Real>0.95912236</Real>
<Real>0.97891223</Real>
<Real>0.97828108</Real>
<Real>0.97720605</Real>
- <Real>0.97623944</Real>
+ <Real>0.97623938</Real>
<Real>0.97527778</Real>
<Real>0.97370166</Real>
<Real>0.97211659</Real>
<Real>0.97065586</Real>
<Real>0.96862853</Real>
<Real>0.96749282</Real>
- <Real>0.96519589</Real>
+ <Real>0.96519583</Real>
<Real>0.96352464</Real>
<Real>0.96173507</Real>
<Real>0.9598251</Real>
<Real>0.94696826</Real>
<Real>0.94605809</Real>
<Real>0.94519788</Real>
- <Real>0.94469965</Real>
+ <Real>0.94469959</Real>
<Real>0.94453776</Real>
<Real>0.94388711</Real>
<Real>0.94351792</Real>
<Real>0.94424397</Real>
<Real>0.94460499</Real>
<Real>0.94499773</Real>
- <Real>0.94563341</Real>
+ <Real>0.94563335</Real>
<Real>0.94637001</Real>
<Real>0.94729853</Real>
<Real>0.94789016</Real>
- <Real>0.9492957</Real>
+ <Real>0.94929564</Real>
<Real>0.95013034</Real>
<Real>0.95128489</Real>
- <Real>0.95289272</Real>
+ <Real>0.95289266</Real>
<Real>0.9538911</Real>
<Real>0.95517474</Real>
<Real>0.95651346</Real>
<Real>0.95739591</Real>
<Real>0.95905453</Real>
- <Real>0.95995909</Real>
- <Real>0.96129471</Real>
+ <Real>0.95995903</Real>
+ <Real>0.96129465</Real>
<Real>0.96225822</Real>
<Real>0.96320486</Real>
<Real>0.96370006</Real>
<Real>0.96682173</Real>
<Real>0.96666759</Real>
<Real>0.96608764</Real>
- <Real>0.96583378</Real>
+ <Real>0.96583372</Real>
<Real>0.96554524</Real>
<Real>0.9645654</Real>
<Real>0.96389264</Real>
<Real>0.94219875</Real>
<Real>0.94211328</Real>
<Real>0.94299912</Real>
- <Real>0.94275576</Real>
+ <Real>0.9427557</Real>
<Real>0.94383633</Real>
<Real>0.94452924</Real>
<Real>0.94524854</Real>
<Real>0.00094952446</Real>
<Real>0.00096428557</Real>
</Sequence>
- <Real Name="Integral">240.60057076136582</Real>
+ <Real Name="Integral">240.60056956927292</Real>
</ReferenceData>
<Real>0.70085573</Real>
<Real>0.63575351</Real>
<Real>0.56621635</Real>
- <Real>0.49021092</Real>
- <Real>0.41346616</Real>
+ <Real>0.49021095</Real>
+ <Real>0.41346619</Real>
<Real>0.32695487</Real>
<Real>0.23947404</Real>
<Real>0.15114491</Real>
- <Real>0.068911292</Real>
- <Real>-0.018343732</Real>
+ <Real>0.068911299</Real>
+ <Real>-0.018343734</Real>
<Real>-0.10222372</Real>
<Real>-0.17851034</Real>
- <Real>-0.2556752</Real>
+ <Real>-0.25567523</Real>
<Real>-0.33233142</Real>
<Real>-0.40906069</Real>
- <Real>-0.47129509</Real>
+ <Real>-0.47129512</Real>
<Real>-0.52546626</Real>
<Real>-0.577564</Real>
<Real>-0.62862134</Real>
<Real>-0.66793829</Real>
<Real>-0.70031446</Real>
<Real>-0.71649045</Real>
- <Real>-0.73871797</Real>
+ <Real>-0.73871803</Real>
<Real>-0.74570811</Real>
<Real>-0.75775838</Real>
<Real>-0.75111836</Real>
<Real>-0.71327382</Real>
<Real>-0.68468153</Real>
<Real>-0.65119308</Real>
- <Real>-0.60643828</Real>
+ <Real>-0.60643834</Real>
<Real>-0.56428039</Real>
<Real>-0.52878588</Real>
<Real>-0.46647781</Real>
<Real>-0.34746397</Real>
<Real>-0.28559166</Real>
<Real>-0.22607948</Real>
- <Real>-0.15521984</Real>
+ <Real>-0.15521985</Real>
<Real>-0.087527938</Real>
- <Real>-0.021420239</Real>
+ <Real>-0.02142024</Real>
<Real>0.042880654</Real>
<Real>0.10628833</Real>
- <Real>0.169659</Real>
+ <Real>0.16965902</Real>
<Real>0.23284768</Real>
- <Real>0.29287171</Real>
+ <Real>0.29287174</Real>
<Real>0.33996987</Real>
<Real>0.3846795</Real>
<Real>0.4326345</Real>
<Real>0.47152123</Real>
<Real>0.50605392</Real>
- <Real>0.53039247</Real>
+ <Real>0.53039253</Real>
<Real>0.552634</Real>
<Real>0.57193202</Real>
<Real>0.57857972</Real>
<Real>0.58289105</Real>
- <Real>0.5867179</Real>
+ <Real>0.58671796</Real>
<Real>0.57858908</Real>
<Real>0.57083881</Real>
<Real>0.54604369</Real>
<Real>0.52723479</Real>
<Real>0.50264382</Real>
<Real>0.4660213</Real>
- <Real>0.42628899</Real>
+ <Real>0.42628902</Real>
<Real>0.38307241</Real>
- <Real>0.33706102</Real>
+ <Real>0.33706105</Real>
<Real>0.29605037</Real>
<Real>0.25652972</Real>
- <Real>0.20322086</Real>
- <Real>0.15179904</Real>
+ <Real>0.20322087</Real>
+ <Real>0.15179905</Real>
<Real>0.1005994</Real>
<Real>0.053715136</Real>
- <Real>-0.0024519849</Real>
+ <Real>-0.0024519851</Real>
<Real>-0.051181242</Real>
- <Real>-0.098006107</Real>
+ <Real>-0.098006114</Real>
<Real>-0.15021549</Real>
- <Real>-0.19438037</Real>
+ <Real>-0.19438039</Real>
<Real>-0.2419212</Real>
- <Real>-0.28681141</Real>
+ <Real>-0.28681144</Real>
<Real>-0.32269096</Real>
<Real>-0.35785615</Real>
- <Real>-0.38878986</Real>
+ <Real>-0.38878989</Real>
<Real>-0.41733515</Real>
- <Real>-0.44138661</Real>
- <Real>-0.44811302</Real>
+ <Real>-0.44138664</Real>
+ <Real>-0.44811305</Real>
<Real>-0.46269119</Real>
<Real>-0.46645734</Real>
- <Real>-0.47144547</Real>
- <Real>-0.47026828</Real>
- <Real>-0.46132591</Real>
+ <Real>-0.4714455</Real>
+ <Real>-0.47026831</Real>
+ <Real>-0.46132594</Real>
<Real>-0.44804862</Real>
<Real>-0.43905067</Real>
<Real>-0.41386977</Real>
- <Real>-0.38871932</Real>
- <Real>-0.36298344</Real>
+ <Real>-0.38871935</Real>
+ <Real>-0.36298347</Real>
<Real>-0.32728449</Real>
<Real>-0.28807318</Real>
<Real>-0.25138158</Real>
- <Real>-0.20809993</Real>
+ <Real>-0.20809995</Real>
<Real>-0.17095725</Real>
<Real>-0.13121751</Real>
<Real>-0.086369112</Real>
<Real>-0.040557165</Real>
<Real>0.0057886061</Real>
- <Real>0.050152354</Real>
- <Real>0.082554683</Real>
- <Real>0.12262377</Real>
+ <Real>0.050152358</Real>
+ <Real>0.082554691</Real>
+ <Real>0.12262378</Real>
<Real>0.15991174</Real>
<Real>0.19901413</Real>
<Real>0.22433081</Real>
<Real>0.37430125</Real>
<Real>0.37234989</Real>
<Real>0.36795545</Real>
- <Real>0.35466814</Real>
- <Real>0.34249428</Real>
+ <Real>0.35466817</Real>
+ <Real>0.34249431</Real>
<Real>0.32283241</Real>
<Real>0.30485228</Real>
<Real>0.28708842</Real>
- <Real>0.25817123</Real>
- <Real>0.22833388</Real>
+ <Real>0.25817126</Real>
+ <Real>0.22833389</Real>
<Real>0.20199309</Real>
- <Real>0.16451561</Real>
- <Real>0.14491151</Real>
+ <Real>0.16451563</Real>
+ <Real>0.14491153</Real>
<Real>0.10192882</Real>
- <Real>0.071294688</Real>
- <Real>0.03886297</Real>
- <Real>0.0038909942</Real>
+ <Real>0.071294695</Real>
+ <Real>0.038862973</Real>
+ <Real>0.0038909945</Real>
<Real>-0.020299936</Real>
<Real>-0.055376183</Real>
<Real>-0.092676423</Real>
<Real>-0.12172776</Real>
<Real>-0.15591128</Real>
- <Real>-0.17435807</Real>
+ <Real>-0.17435808</Real>
<Real>-0.19855703</Real>
<Real>-0.21759604</Real>
- <Real>-0.2342567</Real>
+ <Real>-0.23425671</Real>
<Real>-0.2501609</Real>
- <Real>-0.26639262</Real>
+ <Real>-0.26639265</Real>
<Real>-0.27480033</Real>
<Real>-0.2769354</Real>
<Real>-0.28888217</Real>
<Real>-0.29434624</Real>
- <Real>-0.29110822</Real>
+ <Real>-0.29110825</Real>
<Real>-0.28738263</Real>
<Real>-0.27644756</Real>
<Real>-0.26865044</Real>
- <Real>-0.26017156</Real>
- <Real>-0.24633232</Real>
+ <Real>-0.26017159</Real>
+ <Real>-0.24633233</Real>
<Real>-0.22993366</Real>
- <Real>-0.21036342</Real>
+ <Real>-0.21036343</Real>
<Real>-0.19680807</Real>
<Real>-0.16794714</Real>
- <Real>-0.14982072</Real>
+ <Real>-0.14982073</Real>
<Real>-0.1253354</Real>
- <Real>-0.091874786</Real>
+ <Real>-0.091874793</Real>
<Real>-0.070614971</Real>
<Real>-0.043464452</Real>
<Real>-0.01496042</Real>
<Real>0.0054097059</Real>
- <Real>0.039536409</Real>
+ <Real>0.039536413</Real>
<Real>0.058995754</Real>
<Real>0.087431222</Real>
- <Real>0.10856776</Real>
- <Real>0.12864889</Real>
+ <Real>0.10856777</Real>
+ <Real>0.12864891</Real>
<Real>0.14001396</Real>
<Real>0.15372792</Real>
<Real>0.16890197</Real>
<Real>0.1832682</Real>
- <Real>0.19547278</Real>
+ <Real>0.19547279</Real>
<Real>0.20601526</Real>
- <Real>0.22224233</Real>
+ <Real>0.22224234</Real>
<Real>0.22276327</Real>
- <Real>0.22658545</Real>
+ <Real>0.22658546</Real>
<Real>0.22071175</Real>
<Real>0.21529464</Real>
- <Real>0.21472976</Real>
+ <Real>0.21472977</Real>
<Real>0.20628634</Real>
<Real>0.2043854</Real>
- <Real>0.20099047</Real>
+ <Real>0.20099048</Real>
<Real>0.18452412</Real>
<Real>0.17380331</Real>
<Real>0.15045828</Real>
- <Real>0.13115592</Real>
+ <Real>0.13115594</Real>
<Real>0.1251017</Real>
<Real>0.10887553</Real>
<Real>0.079097189</Real>
<Real>0.061385516</Real>
<Real>0.043582208</Real>
- <Real>0.026112026</Real>
+ <Real>0.026112027</Real>
<Real>0.01028904</Real>
<Real>-0.011082372</Real>
<Real>-0.031485982</Real>
<Real>-0.11735611</Real>
<Real>-0.12887773</Real>
<Real>-0.14123638</Real>
- <Real>-0.1465044</Real>
+ <Real>-0.14650442</Real>
<Real>-0.16701357</Real>
<Real>-0.16918446</Real>
<Real>-0.17741001</Real>
- <Real>-0.18355082</Real>
+ <Real>-0.18355083</Real>
<Real>-0.17439707</Real>
<Real>-0.17841232</Real>
<Real>-0.16846335</Real>
<Real>-0.16889346</Real>
- <Real>-0.17227842</Real>
+ <Real>-0.17227843</Real>
<Real>-0.17070685</Real>
<Real>-0.15092036</Real>
<Real>-0.15334848</Real>
<Real>-0.12954405</Real>
<Real>-0.11295433</Real>
- <Real>-0.09564729</Real>
+ <Real>-0.095647298</Real>
<Real>-0.080270432</Real>
<Real>-0.06479656</Real>
- <Real>-0.055495787</Real>
+ <Real>-0.055495791</Real>
<Real>-0.03828676</Real>
<Real>-0.026680976</Real>
<Real>-0.0092150476</Real>
<Real>-0.0018779703</Real>
- <Real>0.014258126</Real>
+ <Real>0.014258127</Real>
<Real>0.023644734</Real>
- <Real>0.040185757</Real>
+ <Real>0.040185761</Real>
<Real>0.050761189</Real>
<Real>0.060287103</Real>
<Real>0.065608293</Real>
- <Real>0.080362737</Real>
+ <Real>0.080362745</Real>
<Real>0.083356142</Real>
<Real>0.10414003</Real>
<Real>0.11332227</Real>
<Real>0.12000805</Real>
<Real>0.12286837</Real>
- <Real>0.13463314</Real>
- <Real>0.14336015</Real>
+ <Real>0.13463315</Real>
+ <Real>0.14336017</Real>
<Real>0.15017572</Real>
<Real>1.6947948e-05</Real>
<Real>1.720169e-05</Real>
<Real>-0.00016269117</Real>
<Real>-0.00016905692</Real>
</Sequence>
- <Real Name="Integral">0.049695708357404555</Real>
+ <Real Name="Integral">0.049695614293824519</Real>
</ReferenceData>
<Real>0.65127754</Real>
<Real>0.58396322</Real>
<Real>0.50845456</Real>
- <Real>0.43051061</Real>
+ <Real>0.43051064</Real>
<Real>0.35105249</Real>
<Real>0.26587203</Real>
<Real>0.18143819</Real>
<Real>-0.56959647</Real>
<Real>-0.53186774</Real>
<Real>-0.47233477</Real>
- <Real>-0.42289308</Real>
+ <Real>-0.42289311</Real>
<Real>-0.36013916</Real>
<Real>-0.30238953</Real>
<Real>-0.23851043</Real>
<Real>0.033095215</Real>
<Real>0.064681761</Real>
<Real>0.10605592</Real>
- <Real>0.13676096</Real>
+ <Real>0.13676098</Real>
<Real>0.17619975</Real>
<Real>0.20634755</Real>
<Real>0.23424256</Real>
<Real>0.28315401</Real>
<Real>0.31159228</Real>
<Real>0.32296148</Real>
- <Real>0.3337082</Real>
+ <Real>0.33370823</Real>
<Real>0.3472271</Real>
<Real>0.35520938</Real>
<Real>0.35884556</Real>
<Real>-0.24391536</Real>
<Real>-0.23143035</Real>
<Real>-0.21355876</Real>
- <Real>-0.19516303</Real>
+ <Real>-0.19516304</Real>
<Real>-0.17305651</Real>
<Real>-0.14851433</Real>
<Real>-0.13391782</Real>
<Real>-0.0082916291</Real>
<Real>0.028239559</Real>
<Real>0.044842608</Real>
- <Real>0.066723041</Real>
+ <Real>0.066723049</Real>
<Real>0.086818486</Real>
<Real>0.11206019</Real>
<Real>0.13034853</Real>
<Real>0.16507152</Real>
<Real>0.14102207</Real>
<Real>0.13159586</Real>
- <Real>0.11218439</Real>
+ <Real>0.1121844</Real>
<Real>0.090035476</Real>
<Real>0.069503993</Real>
<Real>0.054288179</Real>
<Real>-0.0047988398</Real>
<Real>-0.024740862</Real>
<Real>-0.042112309</Real>
- <Real>-0.057106286</Real>
+ <Real>-0.05710629</Real>
<Real>-0.079133794</Real>
<Real>-0.08531253</Real>
<Real>-0.10183473</Real>
<Real>-8.4526189e-05</Real>
<Real>-0.00019142308</Real>
</Sequence>
- <Real Name="Integral">0.54968381075777017</Real>
+ <Real Name="Integral">0.54968385173596346</Real>
</ReferenceData>
<Real>-0.47069913</Real>
<Real>-0.44109708</Real>
<Real>-0.44801748</Real>
- <Real>-0.44286934</Real>
+ <Real>-0.44286937</Real>
<Real>-0.40633747</Real>
<Real>-0.38159758</Real>
<Real>-0.38471025</Real>
<Real>-0.2954666</Real>
<Real>-0.22040513</Real>
<Real>-0.20538236</Real>
- <Real>-0.17562716</Real>
+ <Real>-0.17562717</Real>
<Real>-0.12331922</Real>
<Real>-0.10807745</Real>
<Real>-0.050659552</Real>
<Real>0.33220929</Real>
<Real>0.33306259</Real>
<Real>0.33927098</Real>
- <Real>0.37198529</Real>
+ <Real>0.37198532</Real>
<Real>0.417128</Real>
<Real>0.4154191</Real>
<Real>0.39847392</Real>
<Real>-0.13059393</Real>
<Real>-0.10386769</Real>
<Real>-0.075353824</Real>
- <Real>-0.065707266</Real>
+ <Real>-0.065707274</Real>
<Real>-0.018131474</Real>
<Real>-0.014584285</Real>
<Real>0.0061852229</Real>
<Real>-0.00067634159</Real>
<Real>-0.0014037765</Real>
</Sequence>
- <Real Name="Integral">-0.13156737015469844</Real>
+ <Real Name="Integral">-0.13156739250644023</Real>
</ReferenceData>
<Real>1</Real>
<Real>1.0026263</Real>
<Real>1.002813</Real>
- <Real>1.0028353</Real>
- <Real>1.00298</Real>
+ <Real>1.0028352</Real>
+ <Real>1.0029799</Real>
<Real>1.0031264</Real>
- <Real>1.0033009</Real>
+ <Real>1.0033008</Real>
<Real>1.0034248</Real>
<Real>1.0034949</Real>
- <Real>1.0036868</Real>
+ <Real>1.0036867</Real>
<Real>1.0038811</Real>
- <Real>1.0039461</Real>
+ <Real>1.0039459</Real>
<Real>1.0039393</Real>
- <Real>1.0039014</Real>
- <Real>1.0041367</Real>
- <Real>1.0040401</Real>
+ <Real>1.0039012</Real>
+ <Real>1.0041366</Real>
+ <Real>1.00404</Real>
<Real>1.0039788</Real>
<Real>1.0041522</Real>
<Real>1.004038</Real>
- <Real>1.00398</Real>
- <Real>1.0038664</Real>
- <Real>1.0038596</Real>
- <Real>1.00374</Real>
- <Real>1.0035481</Real>
- <Real>1.0035183</Real>
+ <Real>1.0039799</Real>
+ <Real>1.0038663</Real>
+ <Real>1.0038595</Real>
+ <Real>1.0037398</Real>
+ <Real>1.003548</Real>
+ <Real>1.0035182</Real>
<Real>1.003296</Real>
<Real>1.0031067</Real>
- <Real>1.0031048</Real>
- <Real>1.0029247</Real>
+ <Real>1.0031047</Real>
+ <Real>1.0029246</Real>
<Real>1.0029281</Real>
- <Real>1.0027844</Real>
+ <Real>1.0027843</Real>
<Real>1.0024828</Real>
<Real>1.0026224</Real>
- <Real>1.0028158</Real>
+ <Real>1.0028157</Real>
<Real>1.0028605</Real>
<Real>1.0029123</Real>
- <Real>1.0031923</Real>
+ <Real>1.0031922</Real>
<Real>1.0032749</Real>
- <Real>1.0033575</Real>
- <Real>1.0034987</Real>
- <Real>1.0035222</Real>
- <Real>1.0038122</Real>
+ <Real>1.0033574</Real>
+ <Real>1.0034986</Real>
+ <Real>1.003522</Real>
+ <Real>1.0038121</Real>
<Real>1.0038934</Real>
<Real>1.0038835</Real>
- <Real>1.003985</Real>
+ <Real>1.0039849</Real>
<Real>1.0040591</Real>
<Real>1.0040226</Real>
<Real>1.004043</Real>
- <Real>1.0039675</Real>
- <Real>1.0040811</Real>
+ <Real>1.0039674</Real>
+ <Real>1.004081</Real>
<Real>1.0039179</Real>
- <Real>1.003891</Real>
+ <Real>1.0038909</Real>
<Real>1.0038245</Real>
<Real>1.0036744</Real>
<Real>1.0038064</Real>
<Real>1.0035872</Real>
<Real>1.003276</Real>
- <Real>1.0033685</Real>
+ <Real>1.0033684</Real>
<Real>1.0032185</Real>
<Real>1.0030403</Real>
<Real>1.0028248</Real>
<Real>1.0028371</Real>
- <Real>1.0028143</Real>
+ <Real>1.0028142</Real>
<Real>1.00263</Real>
- <Real>1.002896</Real>
+ <Real>1.0028958</Real>
<Real>1.0029092</Real>
- <Real>1.0029544</Real>
+ <Real>1.0029542</Real>
<Real>1.0030841</Real>
<Real>1.0033208</Real>
- <Real>1.0033398</Real>
- <Real>1.0034574</Real>
- <Real>1.003602</Real>
- <Real>1.0037988</Real>
+ <Real>1.0033396</Real>
+ <Real>1.0034573</Real>
+ <Real>1.0036019</Real>
+ <Real>1.0037987</Real>
<Real>1.0038297</Real>
<Real>1.0040625</Real>
<Real>1.0039389</Real>
<Real>1.0040179</Real>
<Real>1.0039594</Real>
<Real>1.0039961</Real>
- <Real>1.0039948</Real>
+ <Real>1.0039947</Real>
<Real>1.0040424</Real>
- <Real>1.0039783</Real>
- <Real>1.0039077</Real>
- <Real>1.003881</Real>
- <Real>1.003885</Real>
- <Real>1.0038402</Real>
+ <Real>1.0039781</Real>
+ <Real>1.0039076</Real>
+ <Real>1.0038809</Real>
+ <Real>1.0038849</Real>
+ <Real>1.0038401</Real>
<Real>1.0036588</Real>
<Real>1.0035746</Real>
- <Real>1.0034667</Real>
- <Real>1.0033582</Real>
+ <Real>1.0034666</Real>
+ <Real>1.0033581</Real>
<Real>1.0032606</Real>
- <Real>1.0031298</Real>
- <Real>1.0029531</Real>
- <Real>1.0029762</Real>
+ <Real>1.0031297</Real>
+ <Real>1.0029529</Real>
+ <Real>1.0029761</Real>
<Real>1.0028747</Real>
<Real>1.0028256</Real>
- <Real>1.0029508</Real>
+ <Real>1.0029507</Real>
<Real>1.0031258</Real>
- <Real>1.0032694</Real>
+ <Real>1.0032693</Real>
<Real>1.0033274</Real>
<Real>1.0034236</Real>
<Real>1.0035467</Real>
<Real>1.0037386</Real>
- <Real>1.0037678</Real>
+ <Real>1.0037677</Real>
<Real>1.0038179</Real>
<Real>1.003939</Real>
<Real>1.003979</Real>
<Real>1.00396</Real>
<Real>1.0039968</Real>
- <Real>1.0040765</Real>
+ <Real>1.0040764</Real>
<Real>1.0039867</Real>
<Real>1.0040717</Real>
<Real>1.0039583</Real>
<Real>1.0038927</Real>
- <Real>1.0039638</Real>
+ <Real>1.0039637</Real>
<Real>1.0038794</Real>
<Real>1.003749</Real>
<Real>1.0036877</Real>
- <Real>1.0035897</Real>
- <Real>1.0034268</Real>
+ <Real>1.0035896</Real>
+ <Real>1.0034267</Real>
<Real>1.0034192</Real>
<Real>1.0032343</Real>
- <Real>1.0033185</Real>
- <Real>1.0029762</Real>
+ <Real>1.0033184</Real>
+ <Real>1.0029761</Real>
<Real>1.0031357</Real>
<Real>1.003207</Real>
- <Real>1.0029137</Real>
+ <Real>1.0029136</Real>
<Real>1.0030683</Real>
- <Real>1.0031815</Real>
+ <Real>1.0031813</Real>
<Real>1.0033453</Real>
- <Real>1.0033225</Real>
+ <Real>1.0033224</Real>
<Real>1.0034435</Real>
- <Real>1.0035906</Real>
- <Real>1.0036092</Real>
+ <Real>1.0035905</Real>
+ <Real>1.0036091</Real>
<Real>1.0037911</Real>
- <Real>1.0038091</Real>
- <Real>1.0039287</Real>
+ <Real>1.003809</Real>
+ <Real>1.0039285</Real>
<Real>1.0040594</Real>
<Real>1.003963</Real>
<Real>1.0039438</Real>
- <Real>1.004236</Real>
+ <Real>1.0042359</Real>
<Real>1.0041121</Real>
<Real>1.0039667</Real>
- <Real>1.0039054</Real>
- <Real>1.0039448</Real>
+ <Real>1.0039053</Real>
+ <Real>1.0039446</Real>
<Real>1.0038875</Real>
- <Real>1.0039566</Real>
- <Real>1.0037364</Real>
+ <Real>1.0039564</Real>
+ <Real>1.0037363</Real>
<Real>1.0038453</Real>
- <Real>1.0036628</Real>
+ <Real>1.0036627</Real>
<Real>1.0036311</Real>
- <Real>1.0035908</Real>
+ <Real>1.0035907</Real>
<Real>1.0035497</Real>
- <Real>1.0035697</Real>
- <Real>1.0033984</Real>
+ <Real>1.0035696</Real>
+ <Real>1.0033983</Real>
<Real>1.0034442</Real>
<Real>1.0033016</Real>
- <Real>1.0032938</Real>
- <Real>1.0032138</Real>
+ <Real>1.0032936</Real>
+ <Real>1.0032136</Real>
<Real>1.0033556</Real>
<Real>1.0033462</Real>
<Real>1.0034832</Real>
<Real>1.0035497</Real>
<Real>1.0036097</Real>
<Real>1.003722</Real>
- <Real>1.003818</Real>
+ <Real>1.0038179</Real>
<Real>1.0037833</Real>
- <Real>1.0038013</Real>
- <Real>1.0039788</Real>
+ <Real>1.0038012</Real>
+ <Real>1.0039787</Real>
<Real>1.0039507</Real>
<Real>1.0040078</Real>
- <Real>1.0040756</Real>
- <Real>1.0040264</Real>
+ <Real>1.0040755</Real>
+ <Real>1.0040263</Real>
<Real>1.0038471</Real>
<Real>1.0039368</Real>
<Real>1.0040127</Real>
- <Real>1.0038276</Real>
- <Real>1.0038612</Real>
- <Real>1.0038606</Real>
+ <Real>1.0038275</Real>
+ <Real>1.0038611</Real>
+ <Real>1.0038605</Real>
<Real>1.0038905</Real>
<Real>1.0039259</Real>
- <Real>1.0036864</Real>
+ <Real>1.0036863</Real>
<Real>1.0036433</Real>
- <Real>1.0036565</Real>
+ <Real>1.0036564</Real>
<Real>1.0036999</Real>
<Real>1.003649</Real>
<Real>1.0035149</Real>
<Real>1.003474</Real>
<Real>1.0034895</Real>
<Real>1.003518</Real>
- <Real>1.0036062</Real>
- <Real>1.0036708</Real>
+ <Real>1.0036061</Real>
+ <Real>1.0036707</Real>
<Real>1.0038066</Real>
- <Real>1.0037796</Real>
- <Real>1.0038444</Real>
+ <Real>1.0037795</Real>
+ <Real>1.0038443</Real>
<Real>1.0038437</Real>
<Real>1.0038533</Real>
<Real>1.0040073</Real>
<Real>1.0039427</Real>
- <Real>1.0040594</Real>
+ <Real>1.0040593</Real>
<Real>1.0040133</Real>
<Real>1.0039896</Real>
<Real>1.0041509</Real>
- <Real>1.0040509</Real>
- <Real>1.0040835</Real>
+ <Real>1.0040507</Real>
+ <Real>1.0040834</Real>
<Real>1.0039433</Real>
<Real>1.0039626</Real>
<Real>1.0039417</Real>
<Real>1.0038795</Real>
<Real>1.0038284</Real>
<Real>1.0037864</Real>
- <Real>1.0037988</Real>
- <Real>1.0038092</Real>
- <Real>1.0038675</Real>
- <Real>1.0038058</Real>
+ <Real>1.0037987</Real>
+ <Real>1.0038091</Real>
+ <Real>1.0038674</Real>
+ <Real>1.0038056</Real>
<Real>1.0036768</Real>
- <Real>1.0036308</Real>
+ <Real>1.0036306</Real>
<Real>1.0037299</Real>
- <Real>1.003732</Real>
- <Real>1.0038123</Real>
- <Real>1.0038664</Real>
+ <Real>1.0037318</Real>
+ <Real>1.0038122</Real>
+ <Real>1.0038663</Real>
<Real>1.0037434</Real>
<Real>1.0037843</Real>
- <Real>1.0037801</Real>
+ <Real>1.00378</Real>
<Real>1.0039681</Real>
- <Real>1.003965</Real>
+ <Real>1.0039649</Real>
<Real>1.0039586</Real>
<Real>1.0040132</Real>
<Real>1.0040785</Real>
<Real>1.004027</Real>
<Real>1.0039933</Real>
- <Real>1.0040841</Real>
+ <Real>1.004084</Real>
<Real>1.0040706</Real>
- <Real>1.0040671</Real>
- <Real>1.0040392</Real>
- <Real>1.003906</Real>
- <Real>1.0039387</Real>
- <Real>1.0039802</Real>
+ <Real>1.0040669</Real>
+ <Real>1.004039</Real>
+ <Real>1.0039059</Real>
+ <Real>1.0039386</Real>
+ <Real>1.00398</Real>
<Real>1.0039933</Real>
<Real>1.0039399</Real>
<Real>1.0039351</Real>
- <Real>1.0039452</Real>
+ <Real>1.0039451</Real>
<Real>1.0038861</Real>
- <Real>1.0039293</Real>
+ <Real>1.0039291</Real>
<Real>1.0039737</Real>
- <Real>1.0039887</Real>
+ <Real>1.0039886</Real>
<Real>1.0038258</Real>
<Real>-0.49893174</Real>
<Real>-0.49886957</Real>
<Real>-0.49923372</Real>
<Real>-0.49851903</Real>
</Sequence>
- <Real Name="Integral">127.17098122835159</Real>
+ <Real Name="Integral">127.17096823453903</Real>
</ReferenceData>
<Real>0.65127754</Real>
<Real>0.58396322</Real>
<Real>0.50845456</Real>
- <Real>0.43051061</Real>
+ <Real>0.43051064</Real>
<Real>0.35105249</Real>
<Real>0.26587203</Real>
<Real>0.18143819</Real>
<Real>-0.56959647</Real>
<Real>-0.53186774</Real>
<Real>-0.47233477</Real>
- <Real>-0.42289308</Real>
+ <Real>-0.42289311</Real>
<Real>-0.36013916</Real>
<Real>-0.30238953</Real>
<Real>-0.23851043</Real>
<Real>0.033095215</Real>
<Real>0.064681761</Real>
<Real>0.10605592</Real>
- <Real>0.13676096</Real>
+ <Real>0.13676098</Real>
<Real>0.17619975</Real>
<Real>0.20634755</Real>
<Real>0.23424256</Real>
<Real>0.28315401</Real>
<Real>0.31159228</Real>
<Real>0.32296148</Real>
- <Real>0.3337082</Real>
+ <Real>0.33370823</Real>
<Real>0.3472271</Real>
<Real>0.35520938</Real>
<Real>0.35884556</Real>
<Real>-0.24391536</Real>
<Real>-0.23143035</Real>
<Real>-0.21355876</Real>
- <Real>-0.19516303</Real>
+ <Real>-0.19516304</Real>
<Real>-0.17305651</Real>
<Real>-0.14851433</Real>
<Real>-0.13391782</Real>
<Real>-0.0082916291</Real>
<Real>0.028239559</Real>
<Real>0.044842608</Real>
- <Real>0.066723041</Real>
+ <Real>0.066723049</Real>
<Real>0.086818486</Real>
<Real>0.11206019</Real>
<Real>0.13034853</Real>
<Real>0.16507152</Real>
<Real>0.14102207</Real>
<Real>0.13159586</Real>
- <Real>0.11218439</Real>
+ <Real>0.1121844</Real>
<Real>0.090035476</Real>
<Real>0.069503993</Real>
<Real>0.054288179</Real>
<Real>-0.0047988398</Real>
<Real>-0.024740862</Real>
<Real>-0.042112309</Real>
- <Real>-0.057106286</Real>
+ <Real>-0.05710629</Real>
<Real>-0.079133794</Real>
<Real>-0.08531253</Real>
<Real>-0.10183473</Real>
<Real>-8.4526189e-05</Real>
<Real>-0.00019142308</Real>
</Sequence>
- <Real Name="Integral">0.54968381075777017</Real>
+ <Real Name="Integral">0.54968385173596346</Real>
</ReferenceData>
<Real>0.65127754</Real>
<Real>0.58396322</Real>
<Real>0.50845456</Real>
- <Real>0.43051061</Real>
+ <Real>0.43051064</Real>
<Real>0.35105249</Real>
<Real>0.26587203</Real>
<Real>0.18143819</Real>
<Real>-0.56959647</Real>
<Real>-0.53186774</Real>
<Real>-0.47233477</Real>
- <Real>-0.42289308</Real>
+ <Real>-0.42289311</Real>
<Real>-0.36013916</Real>
<Real>-0.30238953</Real>
<Real>-0.23851043</Real>
<Real>0.033095215</Real>
<Real>0.064681761</Real>
<Real>0.10605592</Real>
- <Real>0.13676096</Real>
+ <Real>0.13676098</Real>
<Real>0.17619975</Real>
<Real>0.20634755</Real>
<Real>0.23424256</Real>
<Real>0.28315401</Real>
<Real>0.31159228</Real>
<Real>0.32296148</Real>
- <Real>0.3337082</Real>
+ <Real>0.33370823</Real>
<Real>0.3472271</Real>
<Real>0.35520938</Real>
<Real>0.35884556</Real>
<Real>-0.24391536</Real>
<Real>-0.23143035</Real>
<Real>-0.21355876</Real>
- <Real>-0.19516303</Real>
+ <Real>-0.19516304</Real>
<Real>-0.17305651</Real>
<Real>-0.14851433</Real>
<Real>-0.13391782</Real>
<Real>-0.0082916291</Real>
<Real>0.028239559</Real>
<Real>0.044842608</Real>
- <Real>0.066723041</Real>
+ <Real>0.066723049</Real>
<Real>0.086818486</Real>
<Real>0.11206019</Real>
<Real>0.13034853</Real>
<Real>0.16507152</Real>
<Real>0.14102207</Real>
<Real>0.13159586</Real>
- <Real>0.11218439</Real>
+ <Real>0.1121844</Real>
<Real>0.090035476</Real>
<Real>0.069503993</Real>
<Real>0.054288179</Real>
<Real>-0.0047988398</Real>
<Real>-0.024740862</Real>
<Real>-0.042112309</Real>
- <Real>-0.057106286</Real>
+ <Real>-0.05710629</Real>
<Real>-0.079133794</Real>
<Real>-0.08531253</Real>
<Real>-0.10183473</Real>
<Real>-8.4526189e-05</Real>
<Real>-0.00019142308</Real>
</Sequence>
- <Real Name="Integral">0.54968381075777017</Real>
+ <Real Name="Integral">0.54968385173596346</Real>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">4</Int>
- <Real>7.8450705641969325</Real>
- <Real>2.6181457557146546</Real>
- <Real>-124.26093946666471</Real>
- <Real>10.991394787222672</Real>
+ <Real>93.478249768350366</Real>
+ <Real>106.54191922850268</Real>
+ <Real>180.89848007846916</Real>
+ <Real>4.0862253105513728</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">3</Int>
- <Real>3.7998256591705597</Real>
- <Real>0.82432886141340755</Real>
- <Real>1.5224332184919882</Real>
+ <Real>0.50082180382778319</Real>
+ <Real>1</Real>
+ <Real>103.32999340041971</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">1</Int>
- <Real>28.790295709638542</Real>
+ <Real>28.633791718066412</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">2</Int>
- <Real>38.159752015476776</Real>
- <Real>0.7956683286568571</Real>
+ <Real>36.978317574617037</Real>
+ <Real>0.80802279377014719</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">5</Int>
- <Real>0.4874505859136089</Real>
- <Real>5.593512789216291</Real>
- <Real>0.5859548866294173</Real>
- <Real>50.639858390563305</Real>
- <Real>-0.002429583757534859</Real>
+ <Real>0.59327596683805051</Real>
+ <Real>16.243870118801162</Real>
+ <Real>0.34166095324008267</Real>
+ <Real>96.46496015267067</Real>
+ <Real>0.014656658285601946</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">7</Int>
- <Real>0.5512395081799536</Real>
- <Real>6.5183422896301373</Real>
- <Real>-0.67296121974350642</Real>
- <Real>31.927951035981206</Real>
- <Real>1.186708933085455</Real>
- <Real>42.008654192804514</Real>
- <Real>-0.0015455112108527936</Real>
+ <Real>0.22969109737961602</Real>
+ <Real>1.6938073214026785</Real>
+ <Real>0.58047615485173931</Real>
+ <Real>23.514547265682463</Real>
+ <Real>0.25468630875040837</Real>
+ <Real>124.69588927433887</Real>
+ <Real>0.0088943612374425204</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">9</Int>
- <Real>15.957817686967116</Real>
- <Real>221006.78653735726</Real>
- <Real>114.88630569000553</Real>
- <Real>-33734.142061366547</Real>
- <Real>7.6914607460936715</Real>
- <Real>1876.2303193471271</Real>
- <Real>0.75265207702909753</Real>
- <Real>24.572552704802391</Real>
- <Real>-138.40801944191745</Real>
+ <Real>0.22969086672835151</Real>
+ <Real>1.6938036292509138</Real>
+ <Real>0.44194868847867191</Real>
+ <Real>23.513035628000434</Real>
+ <Real>0.13852713430709759</Real>
+ <Real>23.519216989391605</Real>
+ <Real>0.25468685075723652</Real>
+ <Real>124.69561162954194</Real>
+ <Real>0.008894417740064841</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">3</Int>
- <Real>49.669368440495347</Real>
- <Real>0.58492964998546348</Real>
- <Real>6.4814246273214478</Real>
+ <Real>6.7235795110785457</Real>
+ <Real>0.42329531456542407</Real>
+ <Real>50.286846364014501</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">6</Int>
- <Real>0.2214172263311579</Real>
- <Real>10.066519430470473</Real>
- <Real>3.7707829256907712</Real>
- <Real>0.72109681793857272</Real>
- <Real>0.49886880345093437</Real>
- <Real>1.000154618776814</Real>
+ <Real>0.49146807677176813</Real>
+ <Real>10.122470863238238</Real>
+ <Real>7.3785890028734542</Real>
+ <Real>1.2022517448320804</Real>
+ <Real>0.41325325292628245</Real>
+ <Real>5.8277122335091622</Real>
</Sequence>
</ReferenceData>
<ReferenceData>
<Sequence Name="result">
<Int Name="Length">2</Int>
- <Real>1.3759034165632673</Real>
- <Real>0.20054113811854915</Real>
+ <Real>0.61188112533079864</Real>
+ <Real>0.085941133945086262</Real>
</Sequence>
</ReferenceData>
--- /dev/null
+0 93.486577
+1 92.437952
+2 92.676899
+3 93.287708
+4 92.524602
+5 93.33982
+6 94.322249
+7 94.152165
+8 93.268429
+9 92.985186
+10 92.447793
+11 93.529461
+12 93.811387
+13 92.128049
+14 93.058745
+15 91.834633
+16 94.592703
+17 94.835931
+18 92.57085
+19 93.350492
+20 93.793811
+21 93.064425
+22 92.358802
+23 93.138607
+24 91.978397
+25 92.883269
+26 94.889234
+27 94.240014
+28 94.141442
+29 94.588848
+30 93.267882
+31 92.642117
+32 94.11224
+33 93.724048
+34 93.050412
+35 92.166011
+36 94.357401
+37 92.104766
+38 93.609402
+39 94.193545
+40 92.616803
+41 93.190927
+42 94.279049
+43 94.066138
+44 92.766857
+45 93.945376
+46 93.410486
+47 92.710437
+48 94.479958
+49 93.37127
+50 94.423499
+51 94.57174
+52 93.645116
+53 93.276059
+54 94.777955
+55 92.01432
+56 93.443601
+57 94.161402
+58 92.767022
+59 94.703138
+60 92.866104
+61 92.709905
+62 93.443922
+63 92.749122
+64 92.996164
+65 93.226291
+66 94.187659
+67 92.239095
+68 93.12234
+69 93.177441
+70 92.728596
+71 92.900251
+72 94.832081
+73 92.066326
+74 93.740011
+75 93.962937
+76 93.559431
+77 92.615057
+78 94.143638
+79 94.78544
+80 93.031341
+81 94.033655
+82 92.418896
+83 92.221132
+84 92.18168
+85 93.919633
+86 94.881543
+87 93.927047
+88 94.674836
+89 93.190343
+90 92.947857
+91 93.298692
+92 94.621239
+93 94.459106
+94 92.906101
+95 94.237764
+96 93.645132
+97 93.660485
+98 92.915948
+99 93.317682
+100 92.120727
+101 93.230415
+102 92.674665
+103 92.159524
+104 93.868182
+105 92.996882
+106 94.363097
+107 94.502451
+108 92.358538
+109 93.717387
+110 92.786443
+111 92.186164
+112 92.229524
+113 92.36136
+114 94.888301
+115 94.425912
+116 92.947455
+117 92.82134
+118 94.760495
+119 94.977061
+120 93.748171
+121 93.605504
+122 94.126006
+123 94.857354
+124 94.460468
+125 92.21732
+126 93.804279
+127 93.52855
+128 94.313085
+129 94.035432
+130 92.303449
+131 94.071923
+132 92.482691
+133 95.028053
+134 92.937195
+135 92.615591
+136 93.580808
+137 94.951146
+138 92.656969
+139 93.256817
+140 92.455751
+141 95.206838
+142 92.984802
+143 92.945197
+144 94.699222
+145 94.56755
+146 95.347739
+147 92.459357
+148 93.171382
+149 95.180269
+150 94.758702
+151 93.379194
+152 93.86194
+153 92.852474
+154 93.279239
+155 93.569225
+156 94.28658
+157 94.625002
+158 93.888688
+159 94.990531
+160 94.31779
+161 93.264266
+162 95.189436
+163 95.923123
+164 93.687356
+165 94.856297
+166 94.09723
+167 96.538322
+168 95.588745
+169 94.174152
+170 94.815074
+171 96.021967
+172 95.894722
+173 97.48244
+174 95.421422
+175 95.738872
+176 96.590132
+177 98.483242
+178 98.857992
+179 97.832445
+180 99.176752
+181 102.0786
+182 101.71236
+183 100.84099
+184 102.06052
+185 103.04261
+186 104.13106
+187 102.48869
+188 103.97843
+189 103.45589
+190 103.46268
+191 103.76516
+192 104.77479
+193 104.83857
+194 103.92016
+195 105.16809
+196 104.87
+197 104.73319
+198 104.75198
+199 106.9468
+200 104.60863
+201 104.743
+202 106.03671
+203 104.92755
+204 104.56418
+205 105.0332
+206 106.65869
+207 104.39706
+208 106.46331
+209 104.91818
+210 106.79197
+211 105.90849
+212 105.78589
+213 106.76522
+214 106.13517
+215 106.84385
+216 107.33566
+217 105.98094
+218 105.46912
+219 104.71926
+220 106.89562
+221 106.95923
+222 106.8085
+223 104.86654
+224 105.4405
+225 106.28891
+226 107.53511
+227 106.43512
+228 105.04485
+229 107.42051
+230 104.91529
+231 105.72297
+232 106.4432
+233 104.84922
+234 107.01712
+235 107.34746
+236 107.16143
+237 106.14437
+238 105.83518
+239 105.85725
+240 106.36342
+241 107.85743
+242 107.0486
+243 105.25835
+244 106.05892
+245 105.21233
+246 105.15961
+247 104.99114
+248 104.98744
+249 106.47118
+250 107.77846
+251 105.95349
+252 106.0645
+253 105.00274
+254 106.7123
+255 106.78758
+256 107.71422
+257 106.74081
+258 107.91838
+259 107.14719
+260 107.57178
+261 107.71677
+262 106.37161
+263 107.76647
+264 105.35725
+265 107.10101
+266 106.95543
+267 106.38286
+268 107.87643
+269 106.18711
+270 105.68704
+271 107.61337
+272 107.30174
+273 107.55044
+274 107.66251
+275 107.97707
+276 108.02862
+277 105.89034
+278 105.8245
+279 106.14353
+280 106.56664
+281 105.80571
+282 107.29199
+283 105.60915
+284 106.63587
+285 106.82268
+286 107.21355
+287 106.76602
+288 106.24447
+289 106.10919
+290 105.73869
+291 105.88343
+292 107.37134
+293 105.92277
+294 107.80709
+295 106.73378
+296 105.28703
+297 107.99916
+298 107.20555
+299 106.72416
+300 106.3867
+301 106.73274
+302 106.99679
+303 107.56717
+304 106.47182
+305 107.21287
+306 106.39991
+307 107.90715
+308 105.09339
+309 107.98186
+310 107.65834
+311 107.39825
+312 107.62341
+313 106.49094
+314 105.92954
+315 107.08428
+316 105.46424
+317 106.1531
+318 105.92754
+319 107.00333
+320 107.87421
+321 105.98247
+322 105.65578
+323 106.11922
+324 105.21628
+325 105.15511
+326 107.23247
+327 107.6705
+328 107.93251
+329 105.8164
+330 105.19615
+331 107.22422
+332 106.46666
+333 105.1668
+334 107.92811
+335 106.38702
+336 106.0796
+337 106.50508
+338 106.15072
+339 105.26819
+340 105.40705
+341 107.25244
+342 105.59738
+343 106.33014
+344 106.44472
+345 106.88026
+346 107.05543
+347 105.61685
+348 106.21632
+349 106.35907
+350 107.49758
+351 106.7012
+352 106.10456
+353 106.93137
+354 107.94109
+355 106.36516
+356 105.22822
+357 107.60685
+358 107.94253
+359 107.26072
+360 105.4247
+361 106.41867
+362 107.69953
+363 105.63072
+364 105.9897
+365 106.46632
+366 106.01961
+367 106.09736
+368 107.5018
+369 106.07223
+370 107.5276
+371 107.19568
+372 107.08564
+373 106.88366
+374 105.74583
+375 107.09943
+376 107.08013
+377 107.45954
+378 107.53612
+379 108.112
+380 107.86845
+381 107.9672
+382 105.9697
+383 107.29257
+384 107.10105
+385 105.21125
+386 107.27416
+387 107.42193
+388 106.61502
+389 106.68435
+390 106.92447
+391 106.36125
+392 106.20247
+393 106.93532
+394 106.22898
+395 106.99196
+396 106.53996
+397 107.53919
+398 105.43453
+399 107.28519
+400 105.87991
--- /dev/null
+1 0.3517423
+17 1.0597503
+33 0.9410293
+49 0.8716023
+65 0.8238423
+81 1.0098463
+97 0.5520603
+113 0.4278303
+129 0.9342923
+145 1.0342603
+161 0.8231863
+177 0.7635853
+193 0.9390683
+209 0.8922163
+225 1.1775333
+241 1.2116213
+257 0.6996663
+273 0.2287793
+289 1.3589063
+305 1.6198323
+321 0.8974273
+337 0.9043013
+353 0.4556703
+369 0.6929363
+385 1.2148863
+401 1.4234283
+417 1.0278523
+433 0.6124313
+449 0.7412923
+465 0.3195843
+481 0.8295103
+497 0.7101183
+513 1.0936703
+529 0.4371843
+545 0.5090993
+561 0.9473693
+577 0.9137463
+593 0.6570103
+609 1.0386173
+625 1.3944973
+641 0.4820373
+657 0.4620413
+673 1.3725933
+689 0.8170593
+705 0.8616223
+721 0.4220783
+737 1.1177493
+753 0.7920653
+769 0.6135533
+785 0.6946533
+801 1.4141893
+817 0.9795273
+833 0.5320023
+849 1.1125843
+865 1.2062423
+881 0.9650543
+897 1.6114013
+913 0.5465363
+929 1.5686003
+945 1.3882573
+961 1.0649153
+977 1.3672613
+993 1.1256913
+1009 0.8110243
+1025 1.0385253
+1041 1.4888273
+1057 1.4508483
+1073 0.9516873
+1089 0.9449663
+1105 0.9077883
+1121 0.5280203
+1137 0.8003353
+1153 0.3715183
+1169 1.1737183
+1185 0.3641943
+1201 0.9924893
+1217 1.4245123
+1233 0.7253083
+1249 1.0137533
+1265 0.8326003
+1281 0.8174633
+1297 0.6439253
+1313 0.9496813
+1329 0.4628573
+1345 1.3082933
+1361 0.3596463
+1377 0.8700153
+1393 1.0316053
+1409 1.2754943
+1425 0.8401303
+1441 1.3051733
+1457 1.7206623
+1473 0.9644593
+1489 0.7674533
+1505 0.9751783
+1521 0.5897953
+1537 0.9886143
+1553 1.5970503
+1569 1.3752943
+1585 1.1988803
+1601 1.0081603
+1617 1.5943493
+1633 0.7405973
+1649 0.6524783
+1665 0.7423293
+1681 1.1931963
+1697 1.2285283
+1713 0.7597933
+1729 0.9291963
+1745 1.5914353
+1761 0.4776883
+1777 1.3845643
+1793 1.4108173
+1809 0.8124203
+1825 1.4891713
+1841 0.8480123
+1857 0.7482803
+1873 1.5186353
+1889 0.8382693
+1905 0.3961233
+1921 0.8167013
+1937 1.3302663
+1953 0.5777863
+1969 1.4406173
+1985 0.4306763
+2001 0.6667753
+2017 -0.0174287
+2033 1.0792663
+2049 0.8624543
+2065 0.7621893
+2081 0.7339673
+2097 0.5926183
+2113 1.0041013
+2129 0.5759933
+2145 1.4805423
+2161 0.9656263
+2177 0.7245373
+2193 0.7458013
+2209 0.7898533
+2225 1.5943193
+2241 0.6717653
+2257 1.2236373
+2273 0.5686843
+2289 1.4059033
+2305 0.8901873
+2321 1.1125003
+2337 0.9947173
+2353 0.7398113
+2369 1.4866303
+2385 1.7170683
+2401 1.1582763
+2417 0.4476063
+2433 0.8380023
+2449 0.7195173
+2465 1.4064993
+2481 1.1095623
+2497 8.13e-05
+2513 0.8899893
+2529 0.9137923
+2545 1.0530973
+2561 1.2563603
+2577 0.3981443
+2593 0.8871053
+2609 0.4778103
+2625 1.1359983
+2641 0.5740863
+2657 1.1778153
+2673 1.3080033
+2689 1.0777173
+2705 1.7451983
+2721 1.9104733
+2737 1.4288603
+2753 2.2134143
+2769 0.9296233
+2785 0.9103213
+2801 0.8418393
+2817 1.2955373
+2833 1.2594723
+2849 0.8617293
+2865 0.2954833
+2881 0.9188513
+2897 1.8631033
+2913 1.9094973
+2929 0.7062733
+2945 1.8554883
+2961 0.9238253
+2977 0.3348663
+2993 0.9353073
+3009 1.2028623
+3025 1.9851123
+3041 1.1466263
+3057 1.5713393
+3073 1.4108013
+3089 0.8622563
+3105 0.8434723
+3121 1.3925823
+3137 1.4874463
+3153 1.2751583
+3169 1.1509373
+3185 1.5736133
+3201 1.5143093
+3217 0.7679413
+3233 1.0632833
+3249 0.9187743
+3265 1.5776033
+3281 1.1580243
+3297 0.8602263
+3313 0.8002063
+3329 0.9148603
+3345 1.3152133
+3361 0.4933133
+3377 0.8761343
+3393 0.6942183
+3409 1.2085083
+3425 1.5972033
+3441 1.2005583
+3457 1.2183353
+3473 0.7577103
+3489 0.6112263
+3505 1.2318773
+3521 1.1326723
+3537 1.1506243
+3553 0.5817003
+3569 0.9014713
+3585 0.8777513
+3601 1.4347963
+3617 0.9277923
+3633 1.1432003
+3649 0.7949343
+3665 1.2829793
+3681 0.2914853
+3697 0.9427083
+3713 1.3510713
+3729 0.8358273
+3745 0.8312423
+3761 1.0880323
+3777 1.7421923
+3793 1.1334733
+3809 0.2839173
+3825 0.8868153
+3841 0.5585983
+3857 0.5283023
+3873 0.2076003
+3889 0.8009313
+3905 0.9836093
+3921 0.8874943
+3937 0.5242053
+3953 0.6777083
+3969 0.2986573
+3985 0.9648863
+4001 0.3936813
+4017 1.5192383
+4033 1.1283383
+4049 0.5402493
+4065 0.8644763
+4081 0.9670683
+4097 1.2198303
+4113 0.8354763
+4129 0.8359803
+4145 1.0664183
+4161 0.4950833
+4177 0.6145223
+4193 0.8997543
+4209 0.2690933
+4225 0.9369703
+4241 1.3704423
+4257 0.5915883
+4273 0.7297033
+4289 0.2907223
+4305 1.0001723
+4321 1.1066633
+4337 1.5133253
+4353 1.3272213
+4369 1.2639593
+4385 0.6724293
+4401 0.6174593
+4417 1.4950913
+4433 0.9158143
+4449 1.9446683
+4465 1.1331683
+4481 0.5225033
+4497 0.8273663
+4513 0.3377963
+4529 0.6311003
+4545 1.4785433
+4561 0.8769733
+4577 0.7451523
+4593 1.0181933
+4609 1.1503113
+4625 1.5157673
+4641 1.2449163
+4657 1.5840353
+4673 1.1181763
+4689 0.8027923
+4705 0.2972533
+4721 0.6047333
+4737 1.3964663
+4753 0.8300983
+4769 1.5832113
+4785 1.7844893
+4801 1.4584933
+4817 1.2985273
+4833 1.1074343
+4849 1.6353573
+4865 1.5705303
+4881 1.7136353
+4897 1.3027233
+4913 1.0545013
+4929 1.0732243
+4945 1.9772463
+4961 1.4843033
+4977 1.4208723
+4993 1.8630953
+5009 0.8611723
+5025 0.7367903
+5041 0.5593233
+5057 0.9989823
+5073 0.5387463
+5089 0.6843843
+5105 0.7769823
+5121 1.1148653
+5137 1.7723363
+5153 0.6389133
+5169 0.6624193
+5185 0.5881323
+5201 0.5588963
+5217 -0.1075927
+5233 0.5866443
+5249 0.4331333
+5265 0.7912493
+5281 0.1760223
+5297 0.8659643
+5313 1.1929063
+5329 0.4952133
+5345 0.6245163
+5361 0.6433383
+5377 0.6050383
+5393 1.1802343
+5409 0.9662523
+5425 1.1449483
+5441 0.4728213
+5457 0.7832533
+5473 1.0352293
+5489 0.8087053
+5505 0.6213123
+5521 0.8185243
+5537 0.2622273
+5553 0.8098043
+5569 0.3489813
+5585 0.5309423
+5601 0.5585833
+5617 -0.0650127
+5633 0.8007403
+5649 1.1783113
+5665 0.8767213
+5681 0.7084473
+5697 0.7087373
+5713 1.2293673
+5729 1.9700823
+5745 1.1370283
+5761 0.8622563
+5777 1.5045673
+5793 1.5535933
+5809 1.6912733
+5825 1.4561283
+5841 1.1133853
+5857 1.0418363
+5873 1.4926883
+5889 0.7510653
+5905 1.4871873
+5921 0.6951643
+5937 0.8722043
+5953 1.0236633
+5969 0.4982503
+5985 1.0484053
+6001 1.4949693
+6017 0.7577103
+6033 1.4756743
+6049 0.9504673
+6065 0.6371503
+6081 0.9289373
+6097 1.6262783
+6113 1.5682493
+6129 1.7294363
+6145 0.8229643
+6161 1.4206823
+6177 1.2592053
+6193 0.8906453
+6209 0.9228333
+6225 0.7128493
+6241 1.0956393
+6257 1.3999223
+6273 1.2278183
+6289 1.3250013
+6305 0.7005813
+6321 0.8609823
+6337 0.9442873
+6353 0.7324343
+6369 1.2111333
+6385 0.9367413
+6401 0.5429203
+6417 0.7491193
+6433 0.6651583
+6449 1.0823103
+6465 1.3247193
+6481 0.8797353
+6497 1.0703243
+6513 0.9107563
+6529 1.8092243
+6545 0.7302523
+6561 1.6127293
+6577 1.2735493
+6593 1.0610703
+6609 0.7459463
+6625 1.0303013
+6641 1.5765423
+6657 1.0391893
+6673 1.1801883
+6689 0.3454793
+6705 0.5029193
+6721 0.7490663
+6737 1.1911363
+6753 1.1398823
+6769 0.8100933
+6785 1.4742703
+6801 1.1287503
+6817 1.4290363
+6833 1.2512253
+6849 1.1904803
+6865 1.1241423
+6881 0.9222003
+6897 0.5952273
+6913 0.8656583
+6929 1.1323363
+6945 1.1635483
+6961 1.3734173
+6977 1.2152753
+6993 0.7563443
+7009 0.6301693
+7025 0.3757603
+7041 0.9500623
+7057 0.9807333
+7073 0.6007513
+7089 1.1412473
+7105 0.8297853
+7121 0.8890813
+7137 0.8713803
+7153 1.0812503
+7169 1.3144803
+7185 1.3724333
+7201 1.7222793
+7217 0.8696183
+7233 1.2517673
+7249 1.0543183
+7265 1.0174603
+7281 1.1306963
+7297 0.4270143
+7313 0.9045993
+7329 1.5317583
+7345 1.1389813
+7361 0.6676073
+7377 0.9130753
+7393 0.5875443
+7409 0.9053693
+7425 0.5835163
+7441 1.1831793
+7457 1.4384203
+7473 0.4676183
+7489 0.7656293
+7505 0.4548773
+7521 1.0648013
+7537 1.0835923
+7553 1.3808793
+7569 1.2125673
+7585 1.2858243
+7601 0.7112093
+7617 0.6617253
+7633 0.8534213
+7649 1.2384923
+7665 1.1525693
+7681 1.8521543
+7697 1.2422913
+7713 1.4445463
+7729 1.0739483
+7745 0.8565343
+7761 0.6647313
+7777 1.2838183
+7793 1.6710863
+7809 1.7629593
+7825 1.7022833
+7841 1.7739233
+7857 1.3768663
+7873 1.7267503
+7889 1.3419463
+7905 1.4587673
+7921 0.8886153
+7937 0.6786163
+7953 0.9597973
+7969 1.7429853
+7985 1.1981473
+8001 1.7303743
+++ /dev/null
-0 0
-1 0.933793
-2 1.64413
-3 2.26881
-4 2.69732
-5 3.18609
-6 3.5633
-7 3.77591
-8 4.083
-9 4.27165
-10 4.50143
-11 4.68664
-12 4.84612
-13 4.93256
-14 5.0487
-15 5.21257
-16 5.29916
-17 5.32257
-18 5.43037
-19 5.50601
-20 5.60239
-21 5.62398
-22 5.72622
-23 5.77308
-24 5.7562
-25 5.84888
-26 5.81699
-27 5.86178
-28 5.94235
-29 5.95981
-30 5.98369
-31 5.95674
-32 6.06792
-33 6.07991
-34 6.11024
-35 6.12377
-36 6.06567
-37 6.16848
-38 6.15328
-39 6.16295
-40 6.16469
-41 6.22697
-42 6.1776
-43 6.27148
-44 6.23865
-45 6.2744
-46 6.27211
-47 6.29498
-48 6.31175
-49 6.28439
-50 6.34758
-51 6.26659
-52 6.29551
-53 6.31016
-54 6.31924
-55 6.32594
-56 6.39353
-57 6.37555
-58 6.4125
-59 6.34882
-60 6.34253
-61 6.37401
-62 6.42251
-63 6.44891
-64 6.41029
-65 6.36925
-66 6.45488
-67 6.47461
-68 6.39323
-69 6.41703
-70 6.47664
-71 6.47864
-72 6.40653
-73 6.4444
-74 6.4901
-75 6.49016
-76 6.43924
-77 6.50607
-78 6.48499
-79 6.4465
-80 6.51673
-81 6.44586
-82 6.53225
-83 6.51432
-84 6.52803
-85 6.49069
-86 6.51123
-87 6.54134
-88 6.50065
-89 6.49403
-90 6.54445
-91 6.56578
-92 6.54692
-93 6.58151
-94 6.5848
-95 6.5796
-96 6.53327
-97 6.57402
-98 6.51995
-99 6.59609
-100 6.59545
-101 6.51956
-102 6.54083
-103 6.53558
-104 6.51182
-105 6.5347
-106 6.55378
-107 6.58284
-108 6.56471
-109 6.57831
-110 6.59092
-111 6.62654
-112 6.54307
-113 6.57412
-114 6.60729
-115 6.59664
-116 6.63216
-117 6.59422
-118 6.61469
-119 6.62636
-120 6.57808
-121 6.57356
-122 6.63797
-123 6.58006
-124 6.62347
-125 6.56521
-126 6.63422
-127 6.57666
-128 6.56846
-129 6.57003
-130 6.56552
-131 6.56682
-132 6.58905
-133 6.59197
-134 6.58329
-135 6.57486
-136 6.6287
-137 6.62418
-138 6.57194
-139 6.62585
-140 6.59216
-141 6.66186
-142 6.61182
-143 6.64493
-144 6.60235
-145 6.64946
-146 6.58789
-147 6.636
-148 6.60824
-149 6.59174
-150 6.67204
-151 6.63844
-152 6.64634
-153 6.61669
-154 6.60086
-155 6.68471
-156 6.6902
-157 6.65257
-158 6.61755
-159 6.67101
-160 6.6674
-161 6.63823
-162 6.63133
-163 6.6367
-164 6.66352
-165 6.696
-166 6.60881
-167 6.67253
-168 6.66735
-169 6.66583
-170 6.66259
-171 6.63691
-172 6.6766
-173 6.6913
-174 6.65969
-175 6.61821
-176 6.6796
-177 6.70715
-178 6.67694
-179 6.69653
-180 6.64676
-181 6.65935
-182 6.6781
-183 6.71181
-184 6.70878
-185 6.67739
-186 6.64651
-187 6.66021
-188 6.71047
-189 6.65987
-190 6.6511
-191 6.71258
-192 6.71689
-193 6.65441
-194 6.65054
-195 6.66806
-196 6.70738
-197 6.64787
-198 6.69967
-199 6.6331
-200 6.63513
-201 6.69851
-202 6.64237
-203 6.70702
-204 6.70229
-205 6.66568
-206 6.6818
-207 6.68629
-208 6.71863
-209 6.73129
-210 6.66952
-211 6.69274
-212 6.67642
-213 6.6782
-214 6.65372
-215 6.7269
-216 6.70297
-217 6.64193
-218 6.65786
-219 6.73524
-220 6.70002
-221 6.6906
-222 6.70867
-223 6.7153
-224 6.66604
-225 6.70376
-226 6.70701
-227 6.71955
-228 6.73206
-229 6.71019
-230 6.70016
-231 6.72956
-232 6.64624
-233 6.73931
-234 6.70839
-235 6.67605
-236 6.67944
-237 6.7271
-238 6.72234
-239 6.66538
-240 6.71936
-241 6.66755
-242 6.68937
-243 6.724
-244 6.68827
-245 6.69197
-246 6.66415
-247 6.69317
-248 6.69017
-249 6.70377
-250 6.70761
-251 6.74691
-252 6.7327
-253 6.66325
-254 6.67565
-255 6.68625
-256 6.71185
-257 6.65938
-258 6.66814
-259 6.69835
-260 6.71306
-261 6.68137
-262 6.68514
-263 6.69503
-264 6.70687
-265 6.70667
-266 6.71226
-267 6.74293
-268 6.6743
-269 6.70273
-270 6.68393
-271 6.71818
-272 6.72623
-273 6.73187
-274 6.73242
-275 6.70873
-276 6.70823
-277 6.69239
-278 6.75188
-279 6.71371
-280 6.7502
-281 6.69877
-282 6.6998
-283 6.66478
-284 6.70866
-285 6.73126
-286 6.73144
-287 6.70046
-288 6.70388
-289 6.6847
-290 6.75264
-291 6.73971
-292 6.70705
-293 6.7407
-294 6.67035
-295 6.72357
-296 6.73306
-297 6.73084
-298 6.68237
-299 6.76651
-300 6.72838
-301 6.76579
-302 6.744
-303 6.6963
-304 6.76216
-305 6.66966
-306 6.67184
-307 6.67226
-308 6.67163
-309 6.68679
-310 6.73726
-311 6.68492
-312 6.73491
-313 6.75372
-314 6.68562
-315 6.73127
-316 6.73742
-317 6.74772
-318 6.71254
-319 6.67357
-320 6.72034
-321 6.71205
-322 6.68542
-323 6.76568
-324 6.69664
-325 6.7009
-326 6.71605
-327 6.73142
-328 6.7524
-329 6.70307
-330 6.7459
-331 6.70831
-332 6.71814
-333 6.76659
-334 6.67721
-335 6.70581
-336 6.75936
-337 6.71626
-338 6.77537
-339 6.73916
-340 6.75086
-341 6.69825
-342 6.73784
-343 6.68206
-344 6.73098
-345 6.69878
-346 6.70428
-347 6.76935
-348 6.68665
-349 6.72539
-350 6.70694
-351 6.683
-352 6.76606
-353 6.74774
-354 6.69002
-355 6.77377
-356 6.71343
-357 6.71858
-358 6.77459
-359 6.72939
-360 6.74046
-361 6.71403
-362 6.70104
-363 6.7676
-364 6.70533
-365 6.68404
-366 6.73426
-367 6.74651
-368 6.73703
-369 6.7093
-370 6.68535
-371 6.78054
-372 6.68628
-373 6.74885
-374 6.75663
-375 6.78102
-376 6.78231
-377 6.69437
-378 6.69903
-379 6.70758
-380 6.74678
-381 6.7422
-382 6.72993
-383 6.68776
-384 6.74261
-385 6.75407
-386 6.76879
-387 6.70534
-388 6.6894
-389 6.71274
-390 6.73197
-391 6.75648
-392 6.74171
-393 6.73153
-394 6.68943
-395 6.75898
-396 6.72838
-397 6.74778
-398 6.72736
-399 6.77489
-400 6.77639
-401 6.70282
-402 6.71161
-403 6.76628
-404 6.77065
-405 6.77603
-406 6.68948
-407 6.76621
-408 6.72207
-409 6.72234
-410 6.75854
-411 6.77273
-412 6.76151
-413 6.74252
-414 6.73153
-415 6.69708
-416 6.74061
-417 6.75818
-418 6.72091
-419 6.74424
-420 6.78673
-421 6.74145
-422 6.69344
-423 6.72471
-424 6.72518
-425 6.76767
-426 6.71893
-427 6.78251
-428 6.78466
-429 6.78094
-430 6.77212
-431 6.74085
-432 6.74283
-433 6.73333
-434 6.78714
-435 6.77238
-436 6.73692
-437 6.73205
-438 6.70148
-439 6.77486
-440 6.76451
-441 6.79155
-442 6.73363
-443 6.69988
-444 6.77442
-445 6.78397
-446 6.7856
-447 6.73778
-448 6.70429
-449 6.7378
-450 6.74829
-451 6.69527
-452 6.75622
-453 6.7867
-454 6.76332
-455 6.72525
-456 6.72494
-457 6.78544
-458 6.7868
-459 6.73629
-460 6.75375
-461 6.71361
-462 6.73621
-463 6.74825
-464 6.70532
-465 6.78782
-466 6.73666
-467 6.79325
-468 6.69962
-469 6.71673
-470 6.72258
-471 6.72758
-472 6.70651
-473 6.70794
-474 6.77938
-475 6.71007
-476 6.7414
-477 6.77997
-478 6.71039
-479 6.77047
-480 6.71826
-481 6.73399
-482 6.69929
-483 6.7929
-484 6.77405
-485 6.79742
-486 6.76882
-487 6.74596
-488 6.71524
-489 6.7947
-490 6.75326
-491 6.72172
-492 6.71176
-493 6.71697
-494 6.74455
-495 6.78932
-496 6.70312
-497 6.72368
-498 6.70806
-499 6.70917
-500 6.79967
--- /dev/null
+0 1.09007
+1 0.877207
+2 0.932819
+3 0.810603
+4 0.791896
+5 0.733364
+6 0.616033
+7 0.604703
+8 0.755853
+9 0.636944
+10 0.715246
+11 0.611031
+12 0.595407
+13 0.63075
+14 0.581412
+15 0.480107
+16 0.51343
+17 0.456788
+18 0.44948
+19 0.435668
+20 0.403863
+21 0.489832
+22 0.495195
+23 0.422594
+24 0.480228
+25 0.482537
+26 0.444401
+27 0.384561
+28 0.34734
+29 0.383506
+30 0.301289
+31 0.382879
+32 0.292707
+33 0.349764
+34 0.302061
+35 0.407986
+36 0.434076
+37 0.335145
+38 0.260494
+39 0.355107
+40 0.25965
+41 0.296108
+42 0.251938
+43 0.255597
+44 0.277503
+45 0.396491
+46 0.28228
+47 0.27551
+48 0.270726
+49 0.292784
+50 0.128415
+51 0.132945
+52 0.23486
+53 0.241426
+54 0.156775
+55 0.2457
+56 0.174162
+57 0.349717
+58 0.186471
+59 0.300814
+60 0.15326
+61 0.231873
+62 0.163605
+63 0.197539
+64 0.260759
+65 0.231805
+66 0.233305
+67 0.135304
+68 0.281354
+69 0.234155
+70 0.176055
+71 0.198226
+72 0.213007
+73 0.225185
+74 0.274865
+75 0.247751
+76 0.105759
+77 0.185042
+78 0.0998593
+79 0.19495
+80 0.0787285
+81 0.31299
+82 0.116831
+83 0.0640667
+84 0.111865
+85 0.127155
+86 0.0863288
+87 0.0788565
+88 0.229011
+89 0.0803561
+90 0.0798449
+91 0.103291
+92 0.148852
+93 0.0735688
+94 0.159655
+95 0.108183
+96 0.0602938
+97 0.148602
+98 0.148881
+99 0.115299
+100 0.0788461
+101 0.181537
+102 0.178337
+103 0.130516
+104 0.0405722
+105 0.178499
+106 0.156721
+107 0.117595
+108 0.190161
+109 0.0353169
+110 0.0107063
+111 0.044126
+112 0.179688
+113 0.0826383
+114 0.161516
+115 0.211201
+116 -0.0253689
+117 0.198065
+118 0.0845166
+119 0.0983947
+120 0.0836944
+121 0.00686212
+122 0.138417
+123 -0.0243976
+124 0.0290124
+125 0.192813
+126 0.190617
+127 0.153298
+128 0.153291
+129 0.160775
+130 0.0652794
+131 0.122779
+132 0.240448
+133 0.072505
+134 0.140685
+135 0.0979626
+136 -0.033676
+137 0.0573126
+138 0.162429
+139 0.0502877
+140 0.176007
+141 0.059104
+142 0.119028
+143 0.0682911
+144 0.0443568
+145 0.0813819
+146 0.046495
+147 0.213421
+148 0.116277
+149 0.211761
+150 0.008448
+151 0.143644
+152 0.203298
+153 0.0675614
+154 0.0783106
+155 0.0764201
+156 0.133036
+157 0.044048
+158 0.0901729
+159 0.00259126
+160 0.088098
+161 0.0468678
+162 0.0648707
+163 0.242762
+164 0.0509743
+165 0.0502586
+166 -0.00215295
+167 0.0977259
+168 0.211509
+169 0.122585
+170 0.115359
+171 0.124501
+172 0.0934846
+173 0.0607657
+174 0.0455892
+175 0.0295404
+176 -0.0230954
+177 0.115156
+178 0.171233
+179 0.0706772
+180 0.0125189
+181 0.0264475
+182 -0.0298884
+183 0.119367
+184 -0.0348317
+185 0.0161417
+186 0.210381
+187 0.15479
+188 0.0625029
+189 0.112174
+190 0.0869022
+191 0.101378
+192 0.0456388
+193 0.0512337
+194 0.0453277
+195 0.141329
+196 0.0583165
+197 -0.100644
+198 -0.0195298
+199 0.0536575
+200 0.0711223
+201 0.0323407
+202 0.0755175
+203 0.151374
+204 0.117336
+205 -0.0224392
+206 0.128654
+207 -0.0316783
+208 -0.00303418
+209 0.15873
+210 -0.0463493
+211 0.027255
+212 0.125538
+213 0.165183
+214 -0.0190088
+215 -0.0792838
+216 0.0241941
+217 0.027236
+218 0.0460868
+219 0.116271
+220 0.00766963
+221 0.0391876
+222 0.011154
+223 0.117667
+224 -0.0332919
+225 0.0515272
+226 -0.0726914
+227 0.0734784
+228 0.00147696
+229 0.133906
+230 0.0234433
+231 0.0647977
+232 0.162547
+233 0.0904404
+234 0.14845
+235 0.00353783
+236 0.175999
+237 0.0489034
+238 0.0733482
+239 0.052908
+240 0.00921323
+241 0.165374
+242 -0.0401146
+243 -0.0183761
+244 0.0339211
+245 -0.0140982
+246 0.162167
+247 0.083416
+248 0.0430574
+249 -0.0532527
+250 0.147883
+251 0.0333986
+252 -0.00200228
+253 0.11477
+254 0.105691
+255 -0.0071386
+256 0.04273
+257 0.0495165
+258 0.0653199
+259 0.0540464
+260 -0.0184714
+261 -0.0314251
+262 -0.0630476
+263 -0.0235271
+264 0.0366917
+265 0.0977608
+266 0.0564583
+267 0.106372
+268 0.0429041
+269 -0.0523068
+270 0.00302942
+271 0.0707971
+272 0.134789
+273 -0.00167497
+274 0.0757873
+275 -0.0220373
+276 0.00421086
+277 0.0529898
+278 0.0759217
+279 0.0753902
+280 0.103891
+281 0.023528
+282 0.0209249
+283 0.117407
+284 -0.00610032
+285 -0.0319453
+286 0.0774823
+287 -0.0464274
+288 0.0475878
+289 -0.0570683
+290 -0.0604167
+291 -0.0104142
+292 0.0735168
+293 0.0105054
+294 0.138414
+295 -0.0477548
+296 -0.0100453
+297 0.0441746
+298 -0.0494035
+299 0.0100155
+300 -0.0115031
+301 -0.0476504
+302 -0.0486697
+303 0.000769082
+304 0.0833754
+305 0.0828655
+306 0.0790107
+307 -0.0217383
+308 -0.0369806
+309 0.0294119
+310 -0.0621266
+311 0.0289071
+312 -0.0395707
+313 0.0415419
+314 0.085237
+315 0.050411
+316 -0.0210736
+317 0.0998225
+318 0.109902
+319 -0.0361716
+320 -0.0101586
+321 -0.0223105
+322 0.00826686
+323 -0.0293934
+324 0.0987037
+325 0.0315682
+326 0.104734
+327 0.0861752
+328 0.0764943
+329 0.000824613
+330 -0.0170506
+331 0.00164694
+332 0.0325245
+333 0.0791064
+334 0.11446
+335 -0.0148082
+336 0.0455031
+337 0.0586395
+338 -0.108668
+339 0.0346456
+340 0.0293549
+341 0.0890408
+342 -0.0106896
+343 0.112917
+344 -0.0296479
+345 0.0522096
+346 0.208064
+347 -0.00224504
+348 0.0121225
+349 -0.105139
+350 -0.00490689
+351 0.0974383
+352 0.045095
+353 0.0802632
+354 -0.000141307
+355 0.00395148
+356 0.156947
+357 0.0607425
+358 0.0246873
+359 0.0152081
+360 0.0806166
+361 0.0525148
+362 0.0390263
+363 -0.082116
+364 -0.0435646
+365 0.0144118
+366 -0.056396
+367 0.0224241
+368 0.0249577
+369 0.032059
+370 0.0710599
+371 0.0617842
+372 0.0870343
+373 -0.0309283
+374 -0.0503928
+375 0.0654705
+376 -0.132297
+377 -0.013171
+378 0.15629
+379 -0.0195602
+380 0.0804519
+381 0.0743904
+382 -0.013602
+383 -0.0495775
+384 -0.0715722
+385 0.150116
+386 0.00559995
+387 0.0252303
+388 0.0791371
+389 -0.00344736
+390 0.051361
+391 -0.0283868
+392 0.00921561
+393 0.0701514
+394 0.110203
+395 -0.00277509
+396 0.155823
+397 0.00737347
+398 -0.00262632
+399 -0.020737
+400 0.0512275
+401 -0.011247
+402 -0.0145163
+403 0.00934317
+404 0.157135
+405 0.102475
+406 0.0686777
+407 0.00617163
+408 -0.024658
+409 0.0290312
+410 0.0888993
+411 0.0475782
+412 -0.0887449
+413 0.068523
+414 -0.029187
+415 0.0134721
+416 0.0268544
+417 0.11216
+418 -0.0280536
+419 -0.0487176
+420 -0.000343852
+421 -0.0204678
+422 0.0346145
+423 -0.0845521
+424 -0.0557614
+425 0.123208
+426 -0.0138119
+427 -0.0674846
+428 -0.0271054
+429 0.0822262
+430 0.0342709
+431 -0.0765876
+432 0.0319721
+433 -0.0739633
+434 0.0233444
+435 -0.00554542
+436 0.0904659
+437 -0.00953448
+438 -0.0605533
+439 -0.0713493
+440 0.017423
+441 0.067253
+442 0.009445
+443 0.00886866
+444 0.0647234
+445 0.0687361
+446 0.00146581
+447 0.0044238
+448 0.0517999
+449 -0.110766
+450 0.0330447
+451 -0.0273471
+452 0.00885882
+453 -0.0352997
+454 0.133558
+455 0.107036
+456 0.0951639
+457 0.0922061
+458 0.107501
+459 0.00342741
+460 0.0689463
+461 -0.109189
+462 0.0251697
+463 -0.104037
+464 0.00550332
+465 0.0253927
+466 0.0473155
+467 -0.0858022
+468 0.0200478
+469 0.0549883
+470 -0.00966798
+471 -0.046954
+472 0.0918711
+473 0.0272196
+474 0.0185346
+475 -0.0329217
+476 -0.0536179
+477 0.0803566
+478 -0.0499015
+479 -0.039256
+480 0.0515844
+481 -0.0254139
+482 0.0432167
+483 0.0460773
+484 -0.0744449
+485 0.134048
+486 -0.0109604
+487 -0.000548119
+488 -0.0285135
+489 0.081197
+490 -0.0220405
+491 0.0489457
+492 0.0368539
+493 0.0850633
+494 0.00580247
+495 0.00387908
+496 -0.046728
+497 0.144958
+498 0.08649
+499 -0.0219127
+500 -0.0263977
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/snprintf.h"
/* must correspond to char *avbar_opt[] declared in main() */
enum {
}
}
-static real anal_ee_inf(double *parm, real T)
+/*! \brief Compute final error estimate.
+ *
+ * See Eqn A17, Hess, JCP 116 (2002) 209-217 for details.
+ */
+static real optimal_error_estimate(double sigma, double fitparm[], real tTotal)
{
- return sqrt(parm[1]*2*parm[0]/T+parm[3]*2*parm[2]/T);
+ double ss = fitparm[1]*fitparm[0]+(1-fitparm[1])*fitparm[2];
+ if ((tTotal <= 0) || (ss <= 0))
+ {
+ fprintf(stderr, "Problem in error estimate: T = %g, ss = %g\n",
+ tTotal, ss);
+ return 0;
+ }
+
+ return sigma*sqrt(2*ss/tTotal);
}
static void estimate_error(const char *eefile, int nb_min, int resol, int n,
double blav, var;
char **leg;
real *tbs, *ybs, rtmp, dens, *fitsig, twooe, tau1_est, tau_sig;
- double fitparm[4];
+ double fitparm[3];
real ee, a, tau1, tau2;
if (n < 4)
*/
fitparm[2] = sqrt(tau1_est*(n-1)*dt);
do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv,
- bDebugMode(), effnERREST, fitparm, 0);
- fitparm[3] = 1-fitparm[1];
+ bDebugMode(), effnERREST, fitparm, 0,
+ NULL);
}
if (bSingleExpFit || fitparm[0] < 0 || fitparm[2] < 0 || fitparm[1] < 0
|| (fitparm[1] > 1 && !bAllowNegLTCorr) || fitparm[2] > (n-1)*dt)
fprintf(stdout, "a fitted parameter is negative\n");
}
fprintf(stdout, "invalid fit: e.e. %g a %g tau1 %g tau2 %g\n",
- sig[s]*anal_ee_inf(fitparm, n*dt),
+ optimal_error_estimate(sig[s], fitparm, n*dt),
fitparm[1], fitparm[0], fitparm[2]);
/* Do a fit with tau2 fixed at the total time.
* One could also choose any other large value for tau2.
fitparm[0] = tau1_est;
fitparm[1] = 0.95;
fitparm[2] = (n-1)*dt;
- fprintf(stderr, "Will fix tau2 at the total time: %g\n", fitparm[2]);
+ fprintf(stdout, "Will fix tau2 at the total time: %g\n", fitparm[2]);
do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv, bDebugMode(),
- effnERREST, fitparm, 4);
- fitparm[3] = 1-fitparm[1];
+ effnERREST, fitparm, 4, NULL);
}
if (bSingleExpFit || fitparm[0] < 0 || fitparm[1] < 0
|| (fitparm[1] > 1 && !bAllowNegLTCorr))
{
fprintf(stdout, "a fitted parameter is negative\n");
fprintf(stdout, "invalid fit: e.e. %g a %g tau1 %g tau2 %g\n",
- sig[s]*anal_ee_inf(fitparm, n*dt),
+ optimal_error_estimate(sig[s], fitparm, n*dt),
fitparm[1], fitparm[0], fitparm[2]);
}
/* Do a single exponential fit */
fitparm[1] = 1.0;
fitparm[2] = 0.0;
do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv, bDebugMode(),
- effnERREST, fitparm, 6);
- fitparm[3] = 1-fitparm[1];
+ effnERREST, fitparm, 6, NULL);
}
}
- ee = sig[s]*anal_ee_inf(fitparm, n*dt);
+ ee = optimal_error_estimate(sig[s], fitparm, n*dt);
a = fitparm[1];
tau1 = fitparm[0];
tau2 = fitparm[2];
if (output_env_get_xvg_format(oenv) == exvgXMGR)
{
fprintf(fp, "@ legend string %d \"av %f\"\n", 2*s, av[s]);
- fprintf(fp, "@ legend string %d \"ee %6g\"\n", 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt));
+ fprintf(fp, "@ legend string %d \"ee %6g\"\n", 2*s+1,
+ optimal_error_estimate(sig[s], fitparm, n*dt));
}
else if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
{
fprintf(fp, "@ s%d legend \"av %f\"\n", 2*s, av[s]);
- fprintf(fp, "@ s%d legend \"ee %6g\"\n", 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt));
+ fprintf(fp, "@ s%d legend \"ee %6g\"\n", 2*s+1,
+ optimal_error_estimate(sig[s], fitparm, n*dt));
}
for (i = 0; i < nbs; i++)
{
ac_fit[1] = 0.95;
ac_fit[2] = 10*acint;
do_lmfit(n/nb_min, ac, fitsig, dt, 0, 0, fitlen*dt, oenv,
- bDebugMode(), effnEXP3, ac_fit, 0);
+ bDebugMode(), effnEXPEXP, ac_fit, 0, NULL);
ac_fit[3] = 1 - ac_fit[1];
fprintf(stdout, "Set %3d: ac erest %g a %g tau1 %g tau2 %g\n",
- s+1, sig[s]*anal_ee_inf(ac_fit, n*dt),
+ s+1, optimal_error_estimate(sig[s], ac_fit, n*dt),
ac_fit[1], ac_fit[0], ac_fit[2]);
fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
static void do_fit(FILE *out, int n, gmx_bool bYdy,
int ny, real *x0, real **val,
- int npargs, t_pargs *ppa, const output_env_t oenv)
+ int npargs, t_pargs *ppa, const output_env_t oenv,
+ const char *fn_fitted)
{
real *c1 = NULL, *sig = NULL;
double *fitparm;
fitparm[0] = 0.5;
fitparm[1] = c1[0];
break;
- case effnEXP3:
+ case effnEXPEXP:
fitparm[0] = 1.0;
fitparm[1] = 0.5*c1[0];
fitparm[2] = 10.0;
fprintf(out, "a%-2d = %12.5e\n", i+1, fitparm[i]);
}
if (do_lmfit(ny, c1, sig, 0, x0, tbeginfit, tendfit,
- oenv, bDebugMode(), efitfn, fitparm, 0))
+ oenv, bDebugMode(), efitfn, fitparm, 0,
+ fn_fitted) > 0)
{
for (i = 0; (i < nparm); i++)
{
xvgrclose(fp);
}
+static void print_fitted_function(const char *fitfile,
+ const char *fn_fitted,
+ gmx_bool bXYdy,
+ int nset,
+ int n,
+ real *t,
+ real **val,
+ int npargs,
+ t_pargs *ppa,
+ output_env_t oenv)
+{
+ FILE *out_fit = gmx_ffopen(fitfile, "w");
+ if (bXYdy && nset >= 2)
+ {
+ do_fit(out_fit, 0, TRUE, n, t, val, npargs, ppa, oenv,
+ fn_fitted);
+ }
+ else
+ {
+ char *buf2 = NULL;
+ int s, buflen = 0;
+ if (NULL != fn_fitted)
+ {
+ buflen = strlen(fn_fitted)+32;
+ snew(buf2, buflen);
+ strncpy(buf2, fn_fitted, buflen);
+ buf2[strlen(buf2)-4] = '\0';
+ }
+ for (s = 0; s < nset; s++)
+ {
+ char *buf = NULL;
+ if (NULL != fn_fitted)
+ {
+ snew(buf, buflen);
+ snprintf(buf, n, "%s_%d.xvg", buf2, s);
+ }
+ do_fit(out_fit, s, FALSE, n, t, val, npargs, ppa, oenv, buf);
+ sfree(buf);
+ }
+ sfree(buf2);
+ }
+ gmx_ffclose(out_fit);
+}
+
int gmx_analyze(int argc, char *argv[])
{
static const char *desc[] = {
"Option [TT]-luzar[tt] performs a Luzar & Chandler kinetics analysis",
"on output from [gmx-hbond]. The input file can be taken directly",
- "from [TT]gmx hbond -ac[tt], and then the same result should be produced."
+ "from [TT]gmx hbond -ac[tt], and then the same result should be produced.[PAR]",
+ "Option [TT]-fitfn[tt] performs curve fitting to a number of different",
+ "curves that make sense in the context of molecular dynamics, mainly",
+ "exponential curves. More information is in the manual. To check the output",
+ "of the fitting procedure the option [TT]-fitted[tt] will print both the",
+ "original data and the fitted function to a new data file. The fitting",
+ "parameters are stored as comment in the output file."
};
static real tb = -1, te = -1, frac = 0.5, filtlen = 0, binwidth = 0.1, aver_start = 0;
static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE;
{ efXVG, "-av", "average", ffOPTWR },
{ efXVG, "-ee", "errest", ffOPTWR },
{ efXVG, "-bal", "ballisitc", ffOPTWR },
+ { efXVG, "-fitted", "fitted", ffOPTWR },
/* { efXVG, "-gem", "geminate", ffOPTWR }, */
{ efLOG, "-g", "fitlog", ffOPTWR }
};
if (fitfile != NULL)
{
- out_fit = gmx_ffopen(fitfile, "w");
- if (bXYdy && nset >= 2)
- {
- do_fit(out_fit, 0, TRUE, n, t, val, npargs, ppa, oenv);
- }
- else
- {
- for (s = 0; s < nset; s++)
- {
- do_fit(out_fit, s, FALSE, n, t, val, npargs, ppa, oenv);
- }
- }
- gmx_ffclose(out_fit);
+ print_fitted_function(fitfile,
+ opt2fn_null("-fitted", NFILE, fnm),
+ bXYdy, nset,
+ n, t, val,
+ npargs, ppa,
+ oenv);
}
printf(" std. dev. relative deviation of\n");
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015, 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.
/*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
/*Fit 1st half of box*/
- do_lmfit(zslices, zDensavg, sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 3);
+ do_lmfit(zslices, zDensavg, sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 8, NULL);
/*Fit 2nd half of box*/
- do_lmfit(zslices, zDensavg, sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 3);
+ do_lmfit(zslices, zDensavg, sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 8, NULL);
/*Initialise the const arrays for storing the average fit parameters*/
avgfit1 = beginfit1;
fit2[k] = avgfit2[k];
}
/*Now fit and store in structures in row-major order int[n][i][j]*/
- do_lmfit(zslices, Densmap[n][i][j], sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, fit1, 1);
+ do_lmfit(zslices, Densmap[n][i][j], sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, fit1, 0, NULL);
int1[n][j+(yslices*i)]->Z = fit1[2];
int1[n][j+(yslices*i)]->t = fit1[3];
- do_lmfit(zslices, Densmap[n][i][j], sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, fit2, 2);
+ do_lmfit(zslices, Densmap[n][i][j], sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, fit2, 0, NULL);
int2[n][j+(yslices*i)]->Z = fit2[2];
int2[n][j+(yslices*i)]->t = fit2[3];
}
integral = print_and_integrate(NULL, calc_nbegin(nx, y[0], tbegin),
dt, y[1], NULL, 1);
integral += do_lmfit(nx, y[1], y[2], dt, y[0], tbegin, tend,
- oenv, TRUE, eFitFn, fitparms, fix);
+ oenv, TRUE, eFitFn, fitparms, fix, NULL);
for (i = 0; i < nx; i++)
{
y[3][i] = fit_function(eFitFn, fitparms, y[0][i]);
}
fitparm[1] = 0.95;
do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(),
- effnEXP2, fitparm, 0);
+ effnEXP2, fitparm, 0, NULL);
}
void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
fitparms[0] = 1;
fitparms[1] = 1;
do_lmfit(ncorr, tcaf[k], sig, dt, 0, 0, ncorr*dt,
- oenv, bDebugMode(), effnVAC, fitparms, 0);
+ oenv, bDebugMode(), effnVAC, fitparms, 0, NULL);
eta = 1000*fitparms[1]*rho/
(4*fitparms[0]*PICO*norm2(kfac[k])/(NANO*NANO));
fprintf(stdout, "k %6.3f tau %6.3f eta %8.5f 10^-3 kg/(m s)\n",
fitparms[0] = 1;
fitparms[1] = 1;
do_lmfit(ncorr, tcafc[k], sig, dt, 0, 0, ncorr*dt,
- oenv, bDebugMode(), effnVAC, fitparms, 0);
+ oenv, bDebugMode(), effnVAC, fitparms, 0, NULL);
eta = 1000*fitparms[1]*rho/
(4*fitparms[0]*PICO*norm2(kfac[kset_c[k]])/(NANO*NANO));
fprintf(stdout,
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff