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