2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/utility/futil.h"
50 #include "gromacs/utility/fatalerror.h"
51 #include "gromacs/math/vec.h"
54 #define MODE(x) ((mode & (x)) == (x))
58 int nrestart
, nout
, P
, fitfn
;
59 gmx_bool bFour
, bNormalize
;
60 real tbeginfit
, tendfit
;
63 static gmx_bool bACFinit
= FALSE
;
70 int sffn2effn(const char **sffn
)
75 for (i
= 0; i
< effnNR
; i
++)
77 if (sffn
[i
+1] && strcmp(sffn
[0], sffn
[i
+1]) == 0)
86 static void low_do_four_core(int nfour
, int nframes
, real c1
[], real cfour
[],
87 int nCos
, gmx_bool bPadding
)
96 for (i
= 0; (i
< nframes
); i
++)
103 for (i
= 0; (i
< nframes
); i
++)
105 cfour
[i
] = cos(c1
[i
]);
109 for (i
= 0; (i
< nframes
); i
++)
111 cfour
[i
] = sin(c1
[i
]);
115 gmx_fatal(FARGS
, "nCos = %d, %s %d", nCos
, __FILE__
, __LINE__
);
117 for (; (i
< nfour
); i
++)
125 /* printf("aver = %g\n",aver); */
126 for (i
= 0; (i
< nframes
); i
++)
133 correl(cfour
-1, cfour
-1, nfour
, ans
-1);
137 for (i
= 0; (i
< nfour
); i
++)
139 cfour
[i
] = ans
[i
]+sqr(aver
);
144 for (i
= 0; (i
< nfour
); i
++)
153 static void do_ac_core(int nframes
, int nout
,
154 real corr
[], real c1
[], int nrestart
,
157 int j
, k
, j3
, jk3
, m
, n
;
163 printf("WARNING: setting number of restarts to 1\n");
169 "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n",
170 nframes
, nout
, nrestart
, mode
);
173 for (j
= 0; (j
< nout
); j
++)
178 /* Loop over starting points. */
179 for (j
= 0; (j
< nframes
); j
+= nrestart
)
183 /* Loop over the correlation length for this starting point */
184 for (k
= 0; (k
< nout
) && (j
+k
< nframes
); k
++)
188 /* Switch over possible ACF types.
189 * It might be more efficient to put the loops inside the switch,
190 * but this is more clear, and save development time!
194 corr
[k
] += c1
[j
]*c1
[j
+k
];
196 else if (MODE(eacCos
))
198 /* Compute the cos (phi(t)-phi(t+dt)) */
199 corr
[k
] += cos(c1
[j
]-c1
[j
+k
]);
201 else if (MODE(eacIden
))
203 /* Check equality (phi(t)==phi(t+dt)) */
204 corr
[k
] += (c1
[j
] == c1
[j
+k
]) ? 1 : 0;
206 else if (MODE(eacP1
) || MODE(eacP2
) || MODE(eacP3
))
208 for (m
= 0; (m
< DIM
); m
++)
213 cth
= cos_angle(xj
, xk
);
215 if (cth
-1.0 > 1.0e-15)
217 printf("j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n",
218 j
, k
, xj
[XX
], xj
[YY
], xj
[ZZ
], xk
[XX
], xk
[YY
], xk
[ZZ
]);
221 corr
[k
] += LegendreP(cth
, mode
); /* 1.5*cth*cth-0.5; */
223 else if (MODE(eacRcross
))
225 for (m
= 0; (m
< DIM
); m
++)
232 corr
[k
] += iprod(rr
, rr
);
234 else if (MODE(eacVector
))
236 for (m
= 0; (m
< DIM
); m
++)
247 gmx_fatal(FARGS
, "\nInvalid mode (%d) in do_ac_core", mode
);
251 /* Correct for the number of points and copy results to the data array */
252 for (j
= 0; (j
< nout
); j
++)
254 n
= (nframes
-j
+(nrestart
-1))/nrestart
;
259 void normalize_acf(int nout
, real corr
[])
266 fprintf(debug
, "Before normalization\n");
267 for (j
= 0; (j
< nout
); j
++)
269 fprintf(debug
, "%5d %10f\n", j
, corr
[j
]);
273 /* Normalisation makes that c[0] = 1.0 and that other points are scaled
284 for (j
= 0; (j
< nout
); j
++)
291 fprintf(debug
, "After normalization\n");
292 for (j
= 0; (j
< nout
); j
++)
294 fprintf(debug
, "%5d %10f\n", j
, corr
[j
]);
299 void average_acf(gmx_bool bVerbose
, int n
, int nitem
, real
**c1
)
306 printf("Averaging correlation functions\n");
309 for (j
= 0; (j
< n
); j
++)
312 for (i
= 0; (i
< nitem
); i
++)
320 void norm_and_scale_vectors(int nframes
, real c1
[], real scale
)
325 for (j
= 0; (j
< nframes
); j
++)
329 for (m
= 0; (m
< DIM
); m
++)
336 void dump_tmp(char *s
, int n
, real c
[])
341 fp
= gmx_ffopen(s
, "w");
342 for (i
= 0; (i
< n
); i
++)
344 fprintf(fp
, "%10d %10g\n", i
, c
[i
]);
349 real
print_and_integrate(FILE *fp
, int n
, real dt
, real c
[], real
*fit
, int nskip
)
354 /* Use trapezoidal rule for calculating integral */
356 for (j
= 0; (j
< n
); j
++)
359 if (fp
&& (nskip
== 0 || j
% nskip
== 0))
361 fprintf(fp
, "%10.3f %10.5f\n", j
*dt
, c0
);
365 sum
+= dt
*(c0
+c
[j
-1]);
373 for (j
= 0; (j
< n
); j
++)
375 if (nskip
== 0 || j
% nskip
== 0)
377 fprintf(fp
, "%10.3f %10.5f\n", j
*dt
, fit
[j
]);
386 real
evaluate_integral(int n
, real x
[], real y
[], real dy
[], real aver_start
,
389 double sum
, sum_var
, w
;
390 double sum_tail
= 0, sum2_tail
= 0;
391 int j
, nsum_tail
= 0;
393 /* Use trapezoidal rule for calculating integral */
396 gmx_fatal(FARGS
, "Evaluating integral: n = %d (file %s, line %d)",
397 n
, __FILE__
, __LINE__
);
402 for (j
= 0; (j
< n
); j
++)
407 w
+= 0.5*(x
[j
] - x
[j
-1]);
411 w
+= 0.5*(x
[j
+1] - x
[j
]);
416 /* Assume all errors are uncorrelated */
417 sum_var
+= sqr(w
*dy
[j
]);
420 if ((aver_start
> 0) && (x
[j
] >= aver_start
))
423 sum2_tail
+= sqrt(sum_var
);
430 sum
= sum_tail
/nsum_tail
;
431 /* This is a worst case estimate, assuming all stddev's are correlated. */
432 *stddev
= sum2_tail
/nsum_tail
;
436 *stddev
= sqrt(sum_var
);
442 void do_four_core(unsigned long mode
, int nfour
, int nf2
, int nframes
,
443 real c1
[], real csum
[], real ctmp
[])
454 /********************************************
456 ********************************************/
457 low_do_four_core(nfour
, nf2
, c1
, csum
, enNorm
, FALSE
);
459 else if (MODE(eacCos
))
461 /***************************************************
463 ***************************************************/
464 /* Copy the data to temp array. Since we need it twice
465 * we can't overwrite original.
467 for (j
= 0; (j
< nf2
); j
++)
472 /* Cosine term of AC function */
473 low_do_four_core(nfour
, nf2
, ctmp
, cfour
, enCos
, FALSE
);
474 for (j
= 0; (j
< nf2
); j
++)
479 /* Sine term of AC function */
480 low_do_four_core(nfour
, nf2
, ctmp
, cfour
, enSin
, FALSE
);
481 for (j
= 0; (j
< nf2
); j
++)
487 else if (MODE(eacP2
))
489 /***************************************************
490 * Legendre polynomials
491 ***************************************************/
492 /* First normalize the vectors */
493 norm_and_scale_vectors(nframes
, c1
, 1.0);
495 /* For P2 thingies we have to do six FFT based correls
496 * First for XX^2, then for YY^2, then for ZZ^2
497 * Then we have to do XY, YZ and XZ (counting these twice)
498 * After that we sum them and normalise
499 * P2(x) = (3 * cos^2 (x) - 1)/2
500 * for unit vectors u and v we compute the cosine as the inner product
501 * cos(u,v) = uX vX + uY vY + uZ vZ
505 * C(t) = | (3 cos^2(u(t'),u(t'+t)) - 1)/2 dt'
510 * P2(u(0),u(t)) = [3 * (uX(0) uX(t) +
512 * uZ(0) uZ(t))^2 - 1]/2
513 * = [3 * ((uX(0) uX(t))^2 +
516 * 2(uX(0) uY(0) uX(t) uY(t)) +
517 * 2(uX(0) uZ(0) uX(t) uZ(t)) +
518 * 2(uY(0) uZ(0) uY(t) uZ(t))) - 1]/2
520 * = [(3/2) * (<uX^2> + <uY^2> + <uZ^2> +
521 * 2<uXuY> + 2<uXuZ> + 2<uYuZ>) - 0.5]
525 /* Because of normalization the number of -0.5 to subtract
526 * depends on the number of data points!
528 for (j
= 0; (j
< nf2
); j
++)
530 csum
[j
] = -0.5*(nf2
-j
);
533 /***** DIAGONAL ELEMENTS ************/
534 for (m
= 0; (m
< DIM
); m
++)
536 /* Copy the vector data in a linear array */
537 for (j
= 0; (j
< nf2
); j
++)
539 ctmp
[j
] = sqr(c1
[DIM
*j
+m
]);
543 sprintf(buf
, "c1diag%d.xvg", m
);
544 dump_tmp(buf
, nf2
, ctmp
);
547 low_do_four_core(nfour
, nf2
, ctmp
, cfour
, enNorm
, FALSE
);
551 sprintf(buf
, "c1dfout%d.xvg", m
);
552 dump_tmp(buf
, nf2
, cfour
);
555 for (j
= 0; (j
< nf2
); j
++)
557 csum
[j
] += fac
*(cfour
[j
]);
560 /******* OFF-DIAGONAL ELEMENTS **********/
561 for (m
= 0; (m
< DIM
); m
++)
563 /* Copy the vector data in a linear array */
565 for (j
= 0; (j
< nf2
); j
++)
567 ctmp
[j
] = c1
[DIM
*j
+m
]*c1
[DIM
*j
+m1
];
572 sprintf(buf
, "c1off%d.xvg", m
);
573 dump_tmp(buf
, nf2
, ctmp
);
575 low_do_four_core(nfour
, nf2
, ctmp
, cfour
, enNorm
, FALSE
);
578 sprintf(buf
, "c1ofout%d.xvg", m
);
579 dump_tmp(buf
, nf2
, cfour
);
582 for (j
= 0; (j
< nf2
); j
++)
584 csum
[j
] += fac
*cfour
[j
];
588 else if (MODE(eacP1
) || MODE(eacVector
))
590 /***************************************************
592 ***************************************************/
595 /* First normalize the vectors */
596 norm_and_scale_vectors(nframes
, c1
, 1.0);
599 /* For vector thingies we have to do three FFT based correls
600 * First for XX, then for YY, then for ZZ
601 * After that we sum them and normalise
603 for (j
= 0; (j
< nf2
); j
++)
607 for (m
= 0; (m
< DIM
); m
++)
609 /* Copy the vector data in a linear array */
610 for (j
= 0; (j
< nf2
); j
++)
612 ctmp
[j
] = c1
[DIM
*j
+m
];
614 low_do_four_core(nfour
, nf2
, ctmp
, cfour
, enNorm
, FALSE
);
615 for (j
= 0; (j
< nf2
); j
++)
623 gmx_fatal(FARGS
, "\nUnknown mode in do_autocorr (%d)", mode
);
627 for (j
= 0; (j
< nf2
); j
++)
629 c1
[j
] = csum
[j
]/(real
)(nframes
-j
);
633 real
fit_acf(int ncorr
, int fitfn
, const output_env_t oenv
, gmx_bool bVerbose
,
634 real tbeginfit
, real tendfit
, real dt
, real c1
[], real
*fit
)
637 real tStart
, tail_corr
, sum
, sumtot
= 0, c_start
, ct_estimate
, *sig
;
638 int i
, j
, jmax
, nf_int
;
641 bPrint
= bVerbose
|| bDebugMode();
652 nf_int
= min(ncorr
, (int)(tendfit
/dt
));
653 sum
= print_and_integrate(debug
, nf_int
, dt
, c1
, NULL
, 1);
655 /* Estimate the correlation time for better fitting */
656 ct_estimate
= 0.5*c1
[0];
657 for (i
= 1; (i
< ncorr
) && (c1
[i
] > 0); i
++)
659 ct_estimate
+= c1
[i
];
661 ct_estimate
*= dt
/c1
[0];
665 printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
666 0.0, dt
*nf_int
, sum
);
667 printf("COR: Relaxation times are computed as fit to an exponential:\n");
668 printf("COR: %s\n", longs_ffn
[fitfn
]);
669 printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit
, min(ncorr
*dt
, tendfit
));
675 printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
676 "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
677 (nfp_ffn
[fitfn
] >= 2) ? " a2 ()" : "",
678 (nfp_ffn
[fitfn
] >= 3) ? " a3 (ps)" : "");
691 for (j
= 0; ((j
< jmax
) && (tStart
< tendfit
) && (tStart
< ncorr
*dt
)); j
++)
693 /* Estimate the correlation time for better fitting */
696 for (i
= 0; (i
< ncorr
) && (dt
*i
< tStart
|| c1
[i
] > 0); i
++)
703 ct_estimate
= 0.5*c1
[i
];
708 ct_estimate
+= c1
[i
];
713 ct_estimate
*= dt
/c_start
;
717 /* The data is strange, so we need to choose somehting */
718 ct_estimate
= tendfit
;
722 fprintf(debug
, "tStart %g ct_estimate: %g\n", tStart
, ct_estimate
);
725 if (fitfn
== effnEXP3
)
727 fitparm
[0] = 0.002*ncorr
*dt
;
729 fitparm
[2] = 0.2*ncorr
*dt
;
733 /* Good initial guess, this increases the probability of convergence */
734 fitparm
[0] = ct_estimate
;
739 /* Generate more or less appropriate sigma's */
740 for (i
= 0; i
< ncorr
; i
++)
742 sig
[i
] = sqrt(ct_estimate
+dt
*i
);
745 nf_int
= min(ncorr
, (int)((tStart
+1e-4)/dt
));
746 sum
= print_and_integrate(debug
, nf_int
, dt
, c1
, NULL
, 1);
747 tail_corr
= do_lmfit(ncorr
, c1
, sig
, dt
, NULL
, tStart
, tendfit
, oenv
,
748 bDebugMode(), fitfn
, fitparm
, 0);
749 sumtot
= sum
+tail_corr
;
750 if (fit
&& ((jmax
== 1) || (j
== 1)))
752 for (i
= 0; (i
< ncorr
); i
++)
754 fit
[i
] = fit_function(fitfn
, fitparm
, i
*dt
);
759 printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart
, sum
, tail_corr
, sumtot
);
760 for (i
= 0; (i
< nfp_ffn
[fitfn
]); i
++)
762 printf(" %11.4e", fitparm
[i
]);
773 void low_do_autocorr(const char *fn
, const output_env_t oenv
, const char *title
,
774 int nframes
, int nitem
, int nout
, real
**c1
,
775 real dt
, unsigned long mode
, int nrestart
,
776 gmx_bool bAver
, gmx_bool bNormalize
,
777 gmx_bool bVerbose
, real tbeginfit
, real tendfit
,
780 FILE *fp
, *gp
= NULL
;
784 real c0
, sum
, Ct2av
, Ctav
;
785 gmx_bool bFour
= acf
.bFour
;
787 /* Check flags and parameters */
788 nout
= get_acfnout();
791 nout
= acf
.nout
= (nframes
+1)/2;
793 else if (nout
> nframes
)
798 if (MODE(eacCos
) && MODE(eacVector
))
800 gmx_fatal(FARGS
, "Incompatible options bCos && bVector (%s, %d)",
803 if ((MODE(eacP3
) || MODE(eacRcross
)) && bFour
)
805 fprintf(stderr
, "Can't combine mode %lu with FFT, turning off FFT\n", mode
);
808 if (MODE(eacNormal
) && MODE(eacVector
))
810 gmx_fatal(FARGS
, "Incompatible mode bits: normal and vector (or Legendre)");
813 /* Print flags and parameters */
816 printf("Will calculate %s of %d thingies for %d frames\n",
817 title
? title
: "autocorrelation", nitem
, nframes
);
818 printf("bAver = %s, bFour = %s bNormalize= %s\n",
819 bool_names
[bAver
], bool_names
[bFour
], bool_names
[bNormalize
]);
820 printf("mode = %lu, dt = %g, nrestart = %d\n", mode
, dt
, nrestart
);
824 c0
= log((double)nframes
)/log(2.0);
834 fprintf(debug
, "Using FFT to calculate %s, #points for FFT = %d\n",
838 /* Allocate temp arrays */
844 nfour
= 0; /* To keep the compiler happy */
849 /* Loop over items (e.g. molecules or dihedrals)
850 * In this loop the actual correlation functions are computed, but without
853 k
= max(1, pow(10, (int)(log(nitem
)/log(100))));
854 for (i
= 0; i
< nitem
; i
++)
856 if (bVerbose
&& ((i
%k
== 0 || i
== nitem
-1)))
858 fprintf(stderr
, "\rThingie %d", i
+1);
863 do_four_core(mode
, nfour
, nframes
, nframes
, c1
[i
], csum
, ctmp
);
867 do_ac_core(nframes
, nout
, ctmp
, c1
[i
], nrestart
, mode
);
872 fprintf(stderr
, "\n");
880 fp
= xvgropen(fn
, title
, "Time (ps)", "C(t)", oenv
);
891 average_acf(bVerbose
, nframes
, nitem
, c1
);
896 normalize_acf(nout
, c1
[0]);
899 if (eFitFn
!= effnNONE
)
901 fit_acf(nout
, eFitFn
, oenv
, fn
!= NULL
, tbeginfit
, tendfit
, dt
, c1
[0], fit
);
902 sum
= print_and_integrate(fp
, nout
, dt
, c1
[0], fit
, 1);
906 sum
= print_and_integrate(fp
, nout
, dt
, c1
[0], NULL
, 1);
909 printf("Correlation time (integral over corrfn): %g (ps)\n", sum
);
915 /* Not averaging. Normalize individual ACFs */
919 gp
= xvgropen("ct-distr.xvg", "Correlation times", "item", "time (ps)", oenv
);
921 for (i
= 0; i
< nitem
; i
++)
925 normalize_acf(nout
, c1
[i
]);
927 if (eFitFn
!= effnNONE
)
929 fit_acf(nout
, eFitFn
, oenv
, fn
!= NULL
, tbeginfit
, tendfit
, dt
, c1
[i
], fit
);
930 sum
= print_and_integrate(fp
, nout
, dt
, c1
[i
], fit
, 1);
934 sum
= print_and_integrate(fp
, nout
, dt
, c1
[i
], NULL
, 1);
938 "CORRelation time (integral over corrfn %d): %g (ps)\n",
946 fprintf(gp
, "%5d %.3f\n", i
, sum
);
957 printf("Average correlation time %.3f Std. Dev. %.3f Error %.3f (ps)\n",
958 Ctav
, sqrt((Ct2av
- sqr(Ctav
))),
959 sqrt((Ct2av
- sqr(Ctav
))/(nitem
-1)));
969 static const char *Leg
[] = { NULL
, "0", "1", "2", "3", NULL
};
971 t_pargs
*add_acf_pargs(int *npargs
, t_pargs
*pa
)
974 { "-acflen", FALSE
, etINT
, {&acf
.nout
},
975 "Length of the ACF, default is half the number of frames" },
976 { "-normalize", FALSE
, etBOOL
, {&acf
.bNormalize
},
978 { "-fftcorr", FALSE
, etBOOL
, {&acf
.bFour
},
979 "HIDDENUse fast fourier transform for correlation function" },
980 { "-nrestart", FALSE
, etINT
, {&acf
.nrestart
},
981 "HIDDENNumber of frames between time origins for ACF when no FFT is used" },
982 { "-P", FALSE
, etENUM
, {Leg
},
983 "Order of Legendre polynomial for ACF (0 indicates none)" },
984 { "-fitfn", FALSE
, etENUM
, {s_ffn
},
986 { "-beginfit", FALSE
, etREAL
, {&acf
.tbeginfit
},
987 "Time where to begin the exponential fit of the correlation function" },
988 { "-endfit", FALSE
, etREAL
, {&acf
.tendfit
},
989 "Time where to end the exponential fit of the correlation function, -1 is until the end" },
991 #define NPA asize(acfpa)
995 snew(ppa
, *npargs
+NPA
);
996 for (i
= 0; (i
< *npargs
); i
++)
1000 for (i
= 0; (i
< NPA
); i
++)
1002 ppa
[*npargs
+i
] = acfpa
[i
];
1010 acf
.fitfn
= effnEXP1
;
1012 acf
.bNormalize
= TRUE
;
1013 acf
.tbeginfit
= 0.0;
1021 void do_autocorr(const char *fn
, const output_env_t oenv
, const char *title
,
1022 int nframes
, int nitem
, real
**c1
,
1023 real dt
, unsigned long mode
, gmx_bool bAver
)
1027 printf("ACF data structures have not been initialised. Call add_acf_pargs\n");
1030 /* Handle enumerated types */
1031 sscanf(Leg
[0], "%d", &acf
.P
);
1032 acf
.fitfn
= sffn2effn(s_ffn
);
1037 mode
= mode
| eacP1
;
1040 mode
= mode
| eacP2
;
1043 mode
= mode
| eacP3
;
1049 low_do_autocorr(fn
, oenv
, title
, nframes
, nitem
, acf
.nout
, c1
, dt
, mode
,
1050 acf
.nrestart
, bAver
, acf
.bNormalize
,
1051 bDebugMode(), acf
.tbeginfit
, acf
.tendfit
,
1055 int get_acfnout(void)
1059 gmx_fatal(FARGS
, "ACF data not initialized yet");
1065 int get_acffitfn(void)
1069 gmx_fatal(FARGS
, "ACF data not initialized yet");
1072 return sffn2effn(s_ffn
);
1075 void cross_corr(int n
, real f
[], real g
[], real corr
[])
1081 correl(f
-1, g
-1, n
, corr
-1);
1085 for (i
= 0; (i
< n
); i
++)
1087 for (j
= i
; (j
< n
); j
++)
1089 corr
[j
-i
] += f
[i
]*g
[j
];
1091 corr
[i
] /= (real
)(n
-i
);