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/basedefinitions.h"
46 #include "gromacs/utility/real.h"
48 static void nrerror(const char error_text[], gmx_bool bExit)
50 fprintf(stderr, "Numerical Recipes run-time error...\n");
51 fprintf(stderr, "%s\n", error_text);
54 fprintf(stderr, "...now exiting to system...\n");
59 /* dont use the keyword vector - it will clash with the
60 * altivec extensions used for powerpc processors.
63 static real *rvector(int nl, int nh)
67 v = (real *)malloc((unsigned) (nh-nl+1)*sizeof(real));
70 nrerror("allocation failure in rvector()", TRUE);
72 /* cppcheck-suppress memleak
73 * free_vector does the same vector arithmetic */
77 static int *ivector(int nl, int nh)
81 v = (int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
84 nrerror("allocation failure in ivector()", TRUE);
86 /* cppcheck-suppress memleak
87 * free_vector does the same vector arithmetic */
92 static real **matrix1(int nrl, int nrh, int ncl, int nch)
97 m = (real **) malloc((unsigned) (nrh-nrl+1)*sizeof(real*));
100 nrerror("allocation failure 1 in matrix1()", TRUE);
104 for (i = nrl; i <= nrh; i++)
106 m[i] = (real *) malloc((unsigned) (nch-ncl+1)*sizeof(real));
109 nrerror("allocation failure 2 in matrix1()", TRUE);
116 static double **dmatrix(int nrl, int nrh, int ncl, int nch)
121 m = (double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
124 nrerror("allocation failure 1 in dmatrix()", TRUE);
128 for (i = nrl; i <= nrh; i++)
130 m[i] = (double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
133 nrerror("allocation failure 2 in dmatrix()", TRUE);
140 static int **imatrix1(int nrl, int nrh, int ncl, int nch)
144 m = (int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
147 nrerror("allocation failure 1 in imatrix1()", TRUE);
151 for (i = nrl; i <= nrh; i++)
153 m[i] = (int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
156 nrerror("allocation failure 2 in imatrix1()", TRUE);
165 static real **submatrix(real **a, int oldrl, int oldrh, int oldcl,
166 int newrl, int newcl)
171 m = (real **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(real*));
174 nrerror("allocation failure in submatrix()", TRUE);
178 for (i = oldrl, j = newrl; i <= oldrh; i++, j++)
180 m[j] = a[i]+oldcl-newcl;
188 static void free_vector(real *v, int nl)
190 free((char*) (v+nl));
193 static void free_ivector(int *v, int nl)
195 free((char*) (v+nl));
198 static void free_dvector(int *v, int nl)
200 free((char*) (v+nl));
205 static void free_matrix(real **m, int nrl, int nrh, int ncl)
209 for (i = nrh; i >= nrl; i--)
211 free((char*) (m[i]+ncl));
213 free((char*) (m+nrl));
216 static real **convert_matrix(real *a, int nrl, int nrh, int ncl, int nch)
218 int i, j, nrow, ncol;
223 m = (real **) malloc((unsigned) (nrow)*sizeof(real*));
226 nrerror("allocation failure in convert_matrix()", TRUE);
229 for (i = 0, j = nrl; i <= nrow-1; i++, j++)
238 static void free_convert_matrix(real **b, int nrl)
240 free((char*) (b+nrl));
243 #define SWAP(a, b) {real temp = (a); (a) = (b); (b) = temp; }
245 static void dump_mat(int n, real **a)
249 for (i = 1; (i <= n); i++)
251 for (j = 1; (j <= n); j++)
253 fprintf(stderr, " %10.3f", a[i][j]);
255 fprintf(stderr, "\n");
259 gmx_bool gaussj(real **a, int n, real **b, int m)
261 int *indxc, *indxr, *ipiv;
262 int i, icol = 0, irow = 0, j, k, l, ll;
263 real big, dum, pivinv;
265 indxc = ivector(1, n);
266 indxr = ivector(1, n);
267 ipiv = ivector(1, n);
268 for (j = 1; j <= n; j++)
272 for (i = 1; i <= n; i++)
275 for (j = 1; j <= n; j++)
279 for (k = 1; k <= n; k++)
283 if (fabs(a[j][k]) >= big)
290 else if (ipiv[k] > 1)
292 nrerror("GAUSSJ: Singular Matrix-1", FALSE);
301 for (l = 1; l <= n; l++)
303 SWAP(a[irow][l], a[icol][l]);
305 for (l = 1; l <= m; l++)
307 SWAP(b[irow][l], b[icol][l]);
312 if (a[icol][icol] == 0.0)
314 fprintf(stderr, "irow = %d, icol = %d\n", irow, icol);
316 nrerror("GAUSSJ: Singular Matrix-2", FALSE);
319 pivinv = 1.0/a[icol][icol];
321 for (l = 1; l <= n; l++)
323 a[icol][l] *= pivinv;
325 for (l = 1; l <= m; l++)
327 b[icol][l] *= pivinv;
329 for (ll = 1; ll <= n; ll++)
335 for (l = 1; l <= n; l++)
337 a[ll][l] -= a[icol][l]*dum;
339 for (l = 1; l <= m; l++)
341 b[ll][l] -= b[icol][l]*dum;
346 for (l = n; l >= 1; l--)
348 if (indxr[l] != indxc[l])
350 for (k = 1; k <= n; k++)
352 SWAP(a[k][indxr[l]], a[k][indxc[l]]);
356 free_ivector(ipiv, 1);
357 free_ivector(indxr, 1);
358 free_ivector(indxc, 1);
366 static void covsrt(real **covar, int ma, int lista[], int mfit)
371 for (j = 1; j < ma; j++)
373 for (i = j+1; i <= ma; i++)
378 for (i = 1; i < mfit; i++)
380 for (j = i+1; j <= mfit; j++)
382 if (lista[j] > lista[i])
384 covar[lista[j]][lista[i]] = covar[i][j];
388 covar[lista[i]][lista[j]] = covar[i][j];
393 for (j = 1; j <= ma; j++)
395 covar[1][j] = covar[j][j];
398 covar[lista[1]][lista[1]] = swap;
399 for (j = 2; j <= mfit; j++)
401 covar[lista[j]][lista[j]] = covar[1][j];
403 for (j = 2; j <= ma; j++)
405 for (i = 1; i <= j-1; i++)
407 covar[i][j] = covar[j][i];
412 #define SWAP(a, b) {swap = (a); (a) = (b); (b) = swap; }
414 static void covsrt_new(real **covar, int ma, int ia[], int mfit)
415 /* Expand in storage the covariance matrix covar, so as to take
416 * into account parameters that are being held fixed. (For the
417 * latter, return zero covariances.)
422 for (i = mfit+1; i <= ma; i++)
424 for (j = 1; j <= i; j++)
426 covar[i][j] = covar[j][i] = 0.0;
430 for (j = ma; j >= 1; j--)
434 for (i = 1; i <= ma; i++)
436 SWAP(covar[i][k], covar[i][j]);
438 for (i = 1; i <= ma; i++)
440 SWAP(covar[k][i], covar[j][i]);
448 static void mrqcof(real x[], real y[], real sig[], int ndata, real a[],
449 int ma, int lista[], int mfit,
450 real **alpha, real beta[], real *chisq,
451 void (*funcs)(real, real *, real *, real *))
454 real ymod, wt, sig2i, dy, *dyda;
456 dyda = rvector(1, ma);
457 for (j = 1; j <= mfit; j++)
459 for (k = 1; k <= j; k++)
466 for (i = 1; i <= ndata; i++)
468 (*funcs)(x[i], a, &ymod, dyda);
469 sig2i = 1.0/(sig[i]*sig[i]);
471 for (j = 1; j <= mfit; j++)
473 wt = dyda[lista[j]]*sig2i;
474 for (k = 1; k <= j; k++)
476 alpha[j][k] += wt*dyda[lista[k]];
480 (*chisq) += dy*dy*sig2i;
482 for (j = 2; j <= mfit; j++)
484 for (k = 1; k <= j-1; k++)
486 alpha[k][j] = alpha[j][k];
489 free_vector(dyda, 1);
493 gmx_bool mrqmin(real x[], real y[], real sig[], int ndata, real a[],
494 int ma, int lista[], int mfit,
495 real **covar, real **alpha, real *chisq,
496 void (*funcs)(real, real *, real *, real *),
500 static real *da, *atry, **oneda, *beta, ochisq;
504 oneda = matrix1(1, mfit, 1, 1);
505 atry = rvector(1, ma);
507 beta = rvector(1, ma);
509 for (j = 1; j <= ma; j++)
512 for (k = 1; k <= mfit; k++)
525 nrerror("Bad LISTA permutation in MRQMIN-1", FALSE);
531 nrerror("Bad LISTA permutation in MRQMIN-2", FALSE);
535 mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, chisq, funcs);
538 for (j = 1; j <= mfit; j++)
540 for (k = 1; k <= mfit; k++)
542 covar[j][k] = alpha[j][k];
544 covar[j][j] = alpha[j][j]*(1.0+(*alamda));
545 oneda[j][1] = beta[j];
547 if (!gaussj(covar, mfit, oneda, 1))
551 for (j = 1; j <= mfit; j++)
557 covsrt(covar, ma, lista, mfit);
558 free_vector(beta, 1);
560 free_vector(atry, 1);
561 free_matrix(oneda, 1, mfit, 1);
564 for (j = 1; j <= ma; j++)
568 for (j = 1; j <= mfit; j++)
570 atry[lista[j]] = a[lista[j]]+da[j];
572 mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, chisq, funcs);
577 for (j = 1; j <= mfit; j++)
579 for (k = 1; k <= mfit; k++)
581 alpha[j][k] = covar[j][k];
584 a[lista[j]] = atry[lista[j]];
596 gmx_bool mrqmin_new(real x[], real y[], real sig[], int ndata, real a[],
597 int ia[], int ma, real **covar, real **alpha, real *chisq,
598 void (*funcs)(real, real [], real *, real []),
600 /* Levenberg-Marquardt method, attempting to reduce the value Chi^2
601 * of a fit between a set of data points x[1..ndata], y[1..ndata]
602 * with individual standard deviations sig[1..ndata], and a nonlinear
603 * function dependent on ma coefficients a[1..ma]. The input array
604 * ia[1..ma] indicates by nonzero entries those components of a that
605 * should be fitted for, and by zero entries those components that
606 * should be held fixed at their input values. The program returns
607 * current best-fit values for the parameters a[1..ma], and
608 * Chi^2 = chisq. The arrays covar[1..ma][1..ma], alpha[1..ma][1..ma]
609 * are used as working space during most iterations. Supply a routine
610 * funcs(x,a,yfit,dyda,ma) that evaluates the fitting function yfit,
611 * and its derivatives dyda[1..ma] with respect to the fitting
612 * parameters a at x. On the first call provide an initial guess for
613 * the parameters a, and set alamda < 0 for initialization (which then
614 * sets alamda=.001). If a step succeeds chisq becomes smaller and
615 * alamda de-creases by a factor of 10. If a step fails alamda grows by
616 * a factor of 10. You must call this routine repeatedly until
617 * convergence is achieved. Then, make one final call with alamda=0,
618 * so that covar[1..ma][1..ma] returns the covariance matrix, and alpha
619 * the curvature matrix.
620 * (Parameters held fixed will return zero covariances.)
623 void covsrt(real **covar, int ma, int ia[], int mfit);
624 gmx_bool gaussj(real **a, int n, real **b, int m);
625 void mrqcof_new(real x[], real y[], real sig[], int ndata, real a[],
626 int ia[], int ma, real **alpha, real beta[], real *chisq,
627 void (*funcs)(real, real [], real *, real []));
630 static real ochisq, *atry, *beta, *da, **oneda;
632 if (*alamda < 0.0) /* Initialization. */
634 atry = rvector(1, ma);
635 beta = rvector(1, ma);
637 for (mfit = 0, j = 1; j <= ma; j++)
644 oneda = matrix1(1, mfit, 1, 1);
646 mrqcof_new(x, y, sig, ndata, a, ia, ma, alpha, beta, chisq, funcs);
648 for (j = 1; j <= ma; j++)
653 for (j = 1; j <= mfit; j++) /* Alter linearized fitting matrix, by augmenting. */
655 for (k = 1; k <= mfit; k++)
657 covar[j][k] = alpha[j][k]; /* diagonal elements. */
659 covar[j][j] = alpha[j][j]*(1.0+(*alamda));
660 oneda[j][1] = beta[j];
662 if (!gaussj(covar, mfit, oneda, 1)) /* Matrix solution. */
666 for (j = 1; j <= mfit; j++)
670 if (*alamda == 0.0) /* Once converged, evaluate covariance matrix. */
672 covsrt_new(covar, ma, ia, mfit);
673 free_matrix(oneda, 1, mfit, 1);
675 free_vector(beta, 1);
676 free_vector(atry, 1);
679 for (j = 0, l = 1; l <= ma; l++) /* Did the trial succeed? */
683 atry[l] = a[l]+da[++j];
686 mrqcof_new(x, y, sig, ndata, atry, ia, ma, covar, da, chisq, funcs);
689 /* Success, accept the new solution. */
692 for (j = 1; j <= mfit; j++)
694 for (k = 1; k <= mfit; k++)
696 alpha[j][k] = covar[j][k];
700 for (l = 1; l <= ma; l++)
705 else /* Failure, increase alamda and return. */
713 void mrqcof_new(real x[], real y[], real sig[], int ndata, real a[],
714 int ia[], int ma, real **alpha, real beta[], real *chisq,
715 void (*funcs)(real, real [], real *, real[]))
716 /* Used by mrqmin to evaluate the linearized fitting matrix alpha, and
717 * vector beta as in (15.5.8), and calculate Chi^2.
720 int i, j, k, l, m, mfit = 0;
721 real ymod, wt, sig2i, dy, *dyda;
723 dyda = rvector(1, ma);
724 for (j = 1; j <= ma; j++)
731 for (j = 1; j <= mfit; j++) /* Initialize (symmetric) alpha), beta. */
733 for (k = 1; k <= j; k++)
740 for (i = 1; i <= ndata; i++) /* Summation loop over all data. */
742 (*funcs)(x[i], a, &ymod, dyda);
743 sig2i = 1.0/(sig[i]*sig[i]);
745 for (j = 0, l = 1; l <= ma; l++)
750 for (j++, k = 0, m = 1; m <= l; m++)
754 alpha[j][++k] += wt*dyda[m];
760 *chisq += dy*dy*sig2i; /* And find Chi^2. */
762 for (j = 2; j <= mfit; j++) /* Fill in the symmetric side. */
764 for (k = 1; k < j; k++)
766 alpha[k][j] = alpha[j][k];
769 free_vector(dyda, 1);