src/gromacs/gmxana/levenmar.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <math.h>
  42 #include <stdio.h>
  43 #include <stdlib.h>
  44
  45 #include "gromacs/utility/basedefinitions.h"
  46 #include "gromacs/utility/real.h"
  47
  48 static void nrerror(const char error_text[], gmx_bool bExit)
  49 {
  50     fprintf(stderr, "Numerical Recipes run-time error...\n");
  51     fprintf(stderr, "%s\n", error_text);
  52     if (bExit)
  53     {
  54         fprintf(stderr, "...now exiting to system...\n");
  55         exit(1);
  56     }
  57 }
  58
  59 /* dont use the keyword vector - it will clash with the
  60  * altivec extensions used for powerpc processors.
  61  */
  62
  63 static real *rvector(int nl, int nh)
  64 {
  65     real *v;
  66
  67     v = (real *)malloc((unsigned) (nh-nl+1)*sizeof(real));
  68     if (!v)
  69     {
  70         nrerror("allocation failure in rvector()", TRUE);
  71     }
  72     /* cppcheck-suppress memleak
  73      * free_vector does the same vector arithmetic */
  74     return v-nl;
  75 }
  76
  77 static int *ivector(int nl, int nh)
  78 {
  79     int *v;
  80
  81     v = (int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
  82     if (!v)
  83     {
  84         nrerror("allocation failure in ivector()", TRUE);
  85     }
  86     /* cppcheck-suppress memleak
  87      * free_vector does the same vector arithmetic */
  88     return v-nl;
  89 }
  90
  91
  92 static real **matrix1(int nrl, int nrh, int ncl, int nch)
  93 {
  94     int    i;
  95     real **m;
  96
  97     m = (real **) malloc((unsigned) (nrh-nrl+1)*sizeof(real*));
  98     if (!m)
  99     {
 100         nrerror("allocation failure 1 in matrix1()", TRUE);
 101     }
 102     m -= nrl;
 103
 104     for (i = nrl; i <= nrh; i++)
 105     {
 106         m[i] = (real *) malloc((unsigned) (nch-ncl+1)*sizeof(real));
 107         if (!m[i])
 108         {
 109             nrerror("allocation failure 2 in matrix1()", TRUE);
 110         }
 111         m[i] -= ncl;
 112     }
 113     return m;
 114 }
 115
 116 static double **dmatrix(int nrl, int nrh, int ncl, int nch)
 117 {
 118     int      i;
 119     double **m;
 120
 121     m = (double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
 122     if (!m)
 123     {
 124         nrerror("allocation failure 1 in dmatrix()", TRUE);
 125     }
 126     m -= nrl;
 127
 128     for (i = nrl; i <= nrh; i++)
 129     {
 130         m[i] = (double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
 131         if (!m[i])
 132         {
 133             nrerror("allocation failure 2 in dmatrix()", TRUE);
 134         }
 135         m[i] -= ncl;
 136     }
 137     return m;
 138 }
 139
 140 static int **imatrix1(int nrl, int nrh, int ncl, int nch)
 141 {
 142     int i, **m;
 143
 144     m = (int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
 145     if (!m)
 146     {
 147         nrerror("allocation failure 1 in imatrix1()", TRUE);
 148     }
 149     m -= nrl;
 150
 151     for (i = nrl; i <= nrh; i++)
 152     {
 153         m[i] = (int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
 154         if (!m[i])
 155         {
 156             nrerror("allocation failure 2 in imatrix1()", TRUE);
 157         }
 158         m[i] -= ncl;
 159     }
 160     return m;
 161 }
 162
 163
 164
 165 static real **submatrix(real **a, int oldrl, int oldrh, int oldcl,
 166                         int newrl, int newcl)
 167 {
 168     int    i, j;
 169     real **m;
 170
 171     m = (real **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(real*));
 172     if (!m)
 173     {
 174         nrerror("allocation failure in submatrix()", TRUE);
 175     }
 176     m -= newrl;
 177
 178     for (i = oldrl, j = newrl; i <= oldrh; i++, j++)
 179     {
 180         m[j] = a[i]+oldcl-newcl;
 181     }
 182
 183     return m;
 184 }
 185
 186
 187
 188 static void free_vector(real *v, int nl)
 189 {
 190     free((char*) (v+nl));
 191 }
 192
 193 static void free_ivector(int *v, int nl)
 194 {
 195     free((char*) (v+nl));
 196 }
 197
 198 static void free_dvector(int *v, int nl)
 199 {
 200     free((char*) (v+nl));
 201 }
 202
 203
 204
 205 static void free_matrix(real **m, int nrl, int nrh, int ncl)
 206 {
 207     int i;
 208
 209     for (i = nrh; i >= nrl; i--)
 210     {
 211         free((char*) (m[i]+ncl));
 212     }
 213     free((char*) (m+nrl));
 214 }
 215
 216 static real **convert_matrix(real *a, int nrl, int nrh, int ncl, int nch)
 217 {
 218     int    i, j, nrow, ncol;
 219     real **m;
 220
 221     nrow = nrh-nrl+1;
 222     ncol = nch-ncl+1;
 223     m    = (real **) malloc((unsigned) (nrow)*sizeof(real*));
 224     if (!m)
 225     {
 226         nrerror("allocation failure in convert_matrix()", TRUE);
 227     }
 228     m -= nrl;
 229     for (i = 0, j = nrl; i <= nrow-1; i++, j++)
 230     {
 231         m[j] = a+ncol*i-ncl;
 232     }
 233     return m;
 234 }
 235
 236
 237
 238 static void free_convert_matrix(real **b, int nrl)
 239 {
 240     free((char*) (b+nrl));
 241 }
 242
 243 #define SWAP(a, b) {real temp = (a); (a) = (b); (b) = temp; }
 244
 245 static void dump_mat(int n, real **a)
 246 {
 247     int i, j;
 248
 249     for (i = 1; (i <= n); i++)
 250     {
 251         for (j = 1; (j <= n); j++)
 252         {
 253             fprintf(stderr, "  %10.3f", a[i][j]);
 254         }
 255         fprintf(stderr, "\n");
 256     }
 257 }
 258
 259 gmx_bool gaussj(real **a, int n, real **b, int m)
 260 {
 261     int *indxc, *indxr, *ipiv;
 262     int  i, icol = 0, irow = 0, j, k, l, ll;
 263     real big, dum, pivinv;
 264
 265     indxc = ivector(1, n);
 266     indxr = ivector(1, n);
 267     ipiv  = ivector(1, n);
 268     for (j = 1; j <= n; j++)
 269     {
 270         ipiv[j] = 0;
 271     }
 272     for (i = 1; i <= n; i++)
 273     {
 274         big = 0.0;
 275         for (j = 1; j <= n; j++)
 276         {
 277             if (ipiv[j] != 1)
 278             {
 279                 for (k = 1; k <= n; k++)
 280                 {
 281                     if (ipiv[k] == 0)
 282                     {
 283                         if (fabs(a[j][k]) >= big)
 284                         {
 285                             big  = fabs(a[j][k]);
 286                             irow = j;
 287                             icol = k;
 288                         }
 289                     }
 290                     else if (ipiv[k] > 1)
 291                     {
 292                         nrerror("GAUSSJ: Singular Matrix-1", FALSE);
 293                         return FALSE;
 294                     }
 295                 }
 296             }
 297         }
 298         ++(ipiv[icol]);
 299         if (irow != icol)
 300         {
 301             for (l = 1; l <= n; l++)
 302             {
 303                 SWAP(a[irow][l], a[icol][l]);
 304             }
 305             for (l = 1; l <= m; l++)
 306             {
 307                 SWAP(b[irow][l], b[icol][l]);
 308             }
 309         }
 310         indxr[i] = irow;
 311         indxc[i] = icol;
 312         if (a[icol][icol] == 0.0)
 313         {
 314             fprintf(stderr, "irow = %d, icol = %d\n", irow, icol);
 315             dump_mat(n, a);
 316             nrerror("GAUSSJ: Singular Matrix-2", FALSE);
 317             return FALSE;
 318         }
 319         pivinv        = 1.0/a[icol][icol];
 320         a[icol][icol] = 1.0;
 321         for (l = 1; l <= n; l++)
 322         {
 323             a[icol][l] *= pivinv;
 324         }
 325         for (l = 1; l <= m; l++)
 326         {
 327             b[icol][l] *= pivinv;
 328         }
 329         for (ll = 1; ll <= n; ll++)
 330         {
 331             if (ll != icol)
 332             {
 333                 dum         = a[ll][icol];
 334                 a[ll][icol] = 0.0;
 335                 for (l = 1; l <= n; l++)
 336                 {
 337                     a[ll][l] -= a[icol][l]*dum;
 338                 }
 339                 for (l = 1; l <= m; l++)
 340                 {
 341                     b[ll][l] -= b[icol][l]*dum;
 342                 }
 343             }
 344         }
 345     }
 346     for (l = n; l >= 1; l--)
 347     {
 348         if (indxr[l] != indxc[l])
 349         {
 350             for (k = 1; k <= n; k++)
 351             {
 352                 SWAP(a[k][indxr[l]], a[k][indxc[l]]);
 353             }
 354         }
 355     }
 356     free_ivector(ipiv, 1);
 357     free_ivector(indxr, 1);
 358     free_ivector(indxc, 1);
 359
 360     return TRUE;
 361 }
 362
 363 #undef SWAP
 364
 365
 366 static void covsrt(real **covar, int ma, int lista[], int mfit)
 367 {
 368     int  i, j;
 369     real swap;
 370
 371     for (j = 1; j < ma; j++)
 372     {
 373         for (i = j+1; i <= ma; i++)
 374         {
 375             covar[i][j] = 0.0;
 376         }
 377     }
 378     for (i = 1; i < mfit; i++)
 379     {
 380         for (j = i+1; j <= mfit; j++)
 381         {
 382             if (lista[j] > lista[i])
 383             {
 384                 covar[lista[j]][lista[i]] = covar[i][j];
 385             }
 386             else
 387             {
 388                 covar[lista[i]][lista[j]] = covar[i][j];
 389             }
 390         }
 391     }
 392     swap = covar[1][1];
 393     for (j = 1; j <= ma; j++)
 394     {
 395         covar[1][j] = covar[j][j];
 396         covar[j][j] = 0.0;
 397     }
 398     covar[lista[1]][lista[1]] = swap;
 399     for (j = 2; j <= mfit; j++)
 400     {
 401         covar[lista[j]][lista[j]] = covar[1][j];
 402     }
 403     for (j = 2; j <= ma; j++)
 404     {
 405         for (i = 1; i <= j-1; i++)
 406         {
 407             covar[i][j] = covar[j][i];
 408         }
 409     }
 410 }
 411
 412 #define SWAP(a, b) {swap = (a); (a) = (b); (b) = swap; }
 413
 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.)
 418  */
 419 {
 420     int  i, j, k;
 421     real swap;
 422     for (i = mfit+1; i <= ma; i++)
 423     {
 424         for (j = 1; j <= i; j++)
 425         {
 426             covar[i][j] = covar[j][i] = 0.0;
 427         }
 428     }
 429     k = mfit;
 430     for (j = ma; j >= 1; j--)
 431     {
 432         if (ia[j])
 433         {
 434             for (i = 1; i <= ma; i++)
 435             {
 436                 SWAP(covar[i][k], covar[i][j]);
 437             }
 438             for (i = 1; i <= ma; i++)
 439             {
 440                 SWAP(covar[k][i], covar[j][i]);
 441             }
 442             k--;
 443         }
 444     }
 445 }
 446 #undef SWAP
 447
 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 *))
 452 {
 453     int  k, j, i;
 454     real ymod, wt, sig2i, dy, *dyda;
 455
 456     dyda = rvector(1, ma);
 457     for (j = 1; j <= mfit; j++)
 458     {
 459         for (k = 1; k <= j; k++)
 460         {
 461             alpha[j][k] = 0.0;
 462         }
 463         beta[j] = 0.0;
 464     }
 465     *chisq = 0.0;
 466     for (i = 1; i <= ndata; i++)
 467     {
 468         (*funcs)(x[i], a, &ymod, dyda);
 469         sig2i = 1.0/(sig[i]*sig[i]);
 470         dy    = y[i]-ymod;
 471         for (j = 1; j <= mfit; j++)
 472         {
 473             wt = dyda[lista[j]]*sig2i;
 474             for (k = 1; k <= j; k++)
 475             {
 476                 alpha[j][k] += wt*dyda[lista[k]];
 477             }
 478             beta[j] += dy*wt;
 479         }
 480         (*chisq) += dy*dy*sig2i;
 481     }
 482     for (j = 2; j <= mfit; j++)
 483     {
 484         for (k = 1; k <= j-1; k++)
 485         {
 486             alpha[k][j] = alpha[j][k];
 487         }
 488     }
 489     free_vector(dyda, 1);
 490 }
 491
 492
 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 *),
 497                 real *alamda)
 498 {
 499     int          k, kk, j, ihit;
 500     static real *da, *atry, **oneda, *beta, ochisq;
 501
 502     if (*alamda < 0.0)
 503     {
 504         oneda = matrix1(1, mfit, 1, 1);
 505         atry  = rvector(1, ma);
 506         da    = rvector(1, ma);
 507         beta  = rvector(1, ma);
 508         kk    = mfit+1;
 509         for (j = 1; j <= ma; j++)
 510         {
 511             ihit = 0;
 512             for (k = 1; k <= mfit; k++)
 513             {
 514                 if (lista[k] == j)
 515                 {
 516                     ihit++;
 517                 }
 518             }
 519             if (ihit == 0)
 520             {
 521                 lista[kk++] = j;
 522             }
 523             else if (ihit > 1)
 524             {
 525                 nrerror("Bad LISTA permutation in MRQMIN-1", FALSE);
 526                 return FALSE;
 527             }
 528         }
 529         if (kk != ma+1)
 530         {
 531             nrerror("Bad LISTA permutation in MRQMIN-2", FALSE);
 532             return FALSE;
 533         }
 534         *alamda = 0.001;
 535         mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, chisq, funcs);
 536         ochisq = (*chisq);
 537     }
 538     for (j = 1; j <= mfit; j++)
 539     {
 540         for (k = 1; k <= mfit; k++)
 541         {
 542             covar[j][k] = alpha[j][k];
 543         }
 544         covar[j][j] = alpha[j][j]*(1.0+(*alamda));
 545         oneda[j][1] = beta[j];
 546     }
 547     if (!gaussj(covar, mfit, oneda, 1))
 548     {
 549         return FALSE;
 550     }
 551     for (j = 1; j <= mfit; j++)
 552     {
 553         da[j] = oneda[j][1];
 554     }
 555     if (*alamda == 0.0)
 556     {
 557         covsrt(covar, ma, lista, mfit);
 558         free_vector(beta, 1);
 559         free_vector(da, 1);
 560         free_vector(atry, 1);
 561         free_matrix(oneda, 1, mfit, 1);
 562         return TRUE;
 563     }
 564     for (j = 1; j <= ma; j++)
 565     {
 566         atry[j] = a[j];
 567     }
 568     for (j = 1; j <= mfit; j++)
 569     {
 570         atry[lista[j]] = a[lista[j]]+da[j];
 571     }
 572     mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, chisq, funcs);
 573     if (*chisq < ochisq)
 574     {
 575         *alamda *= 0.1;
 576         ochisq   = (*chisq);
 577         for (j = 1; j <= mfit; j++)
 578         {
 579             for (k = 1; k <= mfit; k++)
 580             {
 581                 alpha[j][k] = covar[j][k];
 582             }
 583             beta[j]     = da[j];
 584             a[lista[j]] = atry[lista[j]];
 585         }
 586     }
 587     else
 588     {
 589         *alamda *= 10.0;
 590         *chisq   = ochisq;
 591     }
 592     return TRUE;
 593 }
 594
 595
 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 []),
 599                     real *alamda)
 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.)
 621  */
 622 {
 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 []));
 628     int         j, k, l;
 629     static int  mfit;
 630     static real ochisq, *atry, *beta, *da, **oneda;
 631
 632     if (*alamda < 0.0)                    /* Initialization. */
 633     {
 634         atry = rvector(1, ma);
 635         beta = rvector(1, ma);
 636         da   = rvector(1, ma);
 637         for (mfit = 0, j = 1; j <= ma; j++)
 638         {
 639             if (ia[j])
 640             {
 641                 mfit++;
 642             }
 643         }
 644         oneda   = matrix1(1, mfit, 1, 1);
 645         *alamda = 0.001;
 646         mrqcof_new(x, y, sig, ndata, a, ia, ma, alpha, beta, chisq, funcs);
 647         ochisq = (*chisq);
 648         for (j = 1; j <= ma; j++)
 649         {
 650             atry[j] = a[j];
 651         }
 652     }
 653     for (j = 1; j <= mfit; j++) /* Alter linearized fitting matrix, by augmenting. */
 654     {
 655         for (k = 1; k <= mfit; k++)
 656         {
 657             covar[j][k] = alpha[j][k]; /* diagonal elements. */
 658         }
 659         covar[j][j] = alpha[j][j]*(1.0+(*alamda));
 660         oneda[j][1] = beta[j];
 661     }
 662     if (!gaussj(covar, mfit, oneda, 1)) /* Matrix solution. */
 663     {
 664         return FALSE;
 665     }
 666     for (j = 1; j <= mfit; j++)
 667     {
 668         da[j] = oneda[j][1];
 669     }
 670     if (*alamda == 0.0) /* Once converged, evaluate covariance matrix. */
 671     {
 672         covsrt_new(covar, ma, ia, mfit);
 673         free_matrix(oneda, 1, mfit, 1);
 674         free_vector(da, 1);
 675         free_vector(beta, 1);
 676         free_vector(atry, 1);
 677         return TRUE;
 678     }
 679     for (j = 0, l = 1; l <= ma; l++) /* Did the trial succeed? */
 680     {
 681         if (ia[l])
 682         {
 683             atry[l] = a[l]+da[++j];
 684         }
 685     }
 686     mrqcof_new(x, y, sig, ndata, atry, ia, ma, covar, da, chisq, funcs);
 687     if (*chisq < ochisq)
 688     {
 689         /* Success, accept the new solution. */
 690         *alamda *= 0.1;
 691         ochisq   = (*chisq);
 692         for (j = 1; j <= mfit; j++)
 693         {
 694             for (k = 1; k <= mfit; k++)
 695             {
 696                 alpha[j][k] = covar[j][k];
 697             }
 698             beta[j] = da[j];
 699         }
 700         for (l = 1; l <= ma; l++)
 701         {
 702             a[l] = atry[l];
 703         }
 704     }
 705     else                 /* Failure, increase alamda and return. */
 706     {
 707         *alamda *= 10.0;
 708         *chisq   = ochisq;
 709     }
 710     return TRUE;
 711 }
 712
 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.
 718  */
 719 {
 720     int  i, j, k, l, m, mfit = 0;
 721     real ymod, wt, sig2i, dy, *dyda;
 722
 723     dyda = rvector(1, ma);
 724     for (j = 1; j <= ma; j++)
 725     {
 726         if (ia[j])
 727         {
 728             mfit++;
 729         }
 730     }
 731     for (j = 1; j <= mfit; j++) /* Initialize (symmetric) alpha), beta. */
 732     {
 733         for (k = 1; k <= j; k++)
 734         {
 735             alpha[j][k] = 0.0;
 736         }
 737         beta[j] = 0.0;
 738     }
 739     *chisq = 0.0;
 740     for (i = 1; i <= ndata; i++) /* Summation loop over all data. */
 741     {
 742         (*funcs)(x[i], a, &ymod, dyda);
 743         sig2i = 1.0/(sig[i]*sig[i]);
 744         dy    = y[i]-ymod;
 745         for (j = 0, l = 1; l <= ma; l++)
 746         {
 747             if (ia[l])
 748             {
 749                 wt = dyda[l]*sig2i;
 750                 for (j++, k = 0, m = 1; m <= l; m++)
 751                 {
 752                     if (ia[m])
 753                     {
 754                         alpha[j][++k] += wt*dyda[m];
 755                     }
 756                 }
 757                 beta[j] += dy*wt;
 758             }
 759         }
 760         *chisq += dy*dy*sig2i;  /* And find Chi^2. */
 761     }
 762     for (j = 2; j <= mfit; j++) /* Fill in the symmetric side. */
 763     {
 764         for (k = 1; k < j; k++)
 765         {
 766             alpha[k][j] = alpha[j][k];
 767         }
 768     }
 769     free_vector(dyda, 1);
 770 }