Bug Summary

File:gromacs/statistics/statistics.c
Location:line 711, column 67
Description:Array access (from variable 'xr') results in a null pointer dereference

Annotated Source Code

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) 2012,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 "statistics.h"
38#ifdef HAVE_CONFIG_H1
39#include <config.h>
40#endif
41#include <math.h>
42#include "typedefs.h"
43#include "gromacs/utility/smalloc.h"
44#include "gromacs/math/vec.h"
45
46static int gmx_dnint(double x)
47{
48 return (int) (x+0.5);
49}
50
51typedef struct gmx_stats {
52 double aa, a, b, sigma_aa, sigma_a, sigma_b, aver, sigma_aver, error;
53 double rmsd, Rdata, Rfit, Rfitaa, chi2, chi2aa;
54 double *x, *y, *dx, *dy;
55 int computed;
56 int np, np_c, nalloc;
57} gmx_stats;
58
59gmx_stats_t gmx_stats_init()
60{
61 gmx_stats *stats;
62
63 snew(stats, 1)(stats) = save_calloc("stats", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 63, (1), sizeof(*(stats)));
64
65 return (gmx_stats_t) stats;
66}
67
68int gmx_stats_get_npoints(gmx_stats_t gstats, int *N)
69{
70 gmx_stats *stats = (gmx_stats *) gstats;
71
72 *N = stats->np;
73
74 return estatsOK;
75}
76
77int gmx_stats_done(gmx_stats_t gstats)
78{
79 gmx_stats *stats = (gmx_stats *) gstats;
80
81 sfree(stats->x)save_free("stats->x", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 81, (stats->x));
82 stats->x = NULL((void*)0);
83 sfree(stats->y)save_free("stats->y", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 83, (stats->y));
84 stats->y = NULL((void*)0);
85 sfree(stats->dx)save_free("stats->dx", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 85, (stats->dx));
86 stats->dx = NULL((void*)0);
87 sfree(stats->dy)save_free("stats->dy", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 87, (stats->dy));
88 stats->dy = NULL((void*)0);
89
90 return estatsOK;
91}
92
93int gmx_stats_add_point(gmx_stats_t gstats, double x, double y,
94 double dx, double dy)
95{
96 gmx_stats *stats = (gmx_stats *) gstats;
97 int i;
98
99 if (stats->np+1 >= stats->nalloc)
100 {
101 if (stats->nalloc == 0)
102 {
103 stats->nalloc = 1024;
104 }
105 else
106 {
107 stats->nalloc *= 2;
108 }
109 srenew(stats->x, stats->nalloc)(stats->x) = save_realloc("stats->x", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 109, (stats->x), (stats->nalloc), sizeof(*(stats->
x)));
110 srenew(stats->y, stats->nalloc)(stats->y) = save_realloc("stats->y", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 110, (stats->y), (stats->nalloc), sizeof(*(stats->
y)));
111 srenew(stats->dx, stats->nalloc)(stats->dx) = save_realloc("stats->dx", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 111, (stats->dx), (stats->nalloc), sizeof(*(stats->
dx)));
112 srenew(stats->dy, stats->nalloc)(stats->dy) = save_realloc("stats->dy", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 112, (stats->dy), (stats->nalloc), sizeof(*(stats->
dy)));
113 for (i = stats->np; (i < stats->nalloc); i++)
114 {
115 stats->x[i] = 0;
116 stats->y[i] = 0;
117 stats->dx[i] = 0;
118 stats->dy[i] = 0;
119 }
120 }
121 stats->x[stats->np] = x;
122 stats->y[stats->np] = y;
123 stats->dx[stats->np] = dx;
124 stats->dy[stats->np] = dy;
125 stats->np++;
126 stats->computed = 0;
127
128 return estatsOK;
129}
130
131int gmx_stats_get_point(gmx_stats_t gstats, real *x, real *y,
132 real *dx, real *dy, real level)
133{
134 gmx_stats *stats = (gmx_stats *) gstats;
135 int ok, outlier;
136 real rmsd, r;
137
138 if ((ok = gmx_stats_get_rmsd(gstats, &rmsd)) != estatsOK)
139 {
140 return ok;
141 }
142 outlier = 0;
143 while ((outlier == 0) && (stats->np_c < stats->np))
144 {
145 r = fabs(stats->x[stats->np_c] - stats->y[stats->np_c]);
146 outlier = (r > rmsd*level);
147 if (outlier)
148 {
149 if (NULL((void*)0) != x)
150 {
151 *x = stats->x[stats->np_c];
152 }
153 if (NULL((void*)0) != y)
154 {
155 *y = stats->y[stats->np_c];
156 }
157 if (NULL((void*)0) != dx)
158 {
159 *dx = stats->dx[stats->np_c];
160 }
161 if (NULL((void*)0) != dy)
162 {
163 *dy = stats->dy[stats->np_c];
164 }
165 }
166 stats->np_c++;
167
168 if (outlier)
169 {
170 return estatsOK;
171 }
172 }
173
174 stats->np_c = 0;
175
176 return estatsNO_POINTS;
177}
178
179int gmx_stats_add_points(gmx_stats_t gstats, int n, real *x, real *y,
180 real *dx, real *dy)
181{
182 int i, ok;
183
184 for (i = 0; (i < n); i++)
185 {
186 if ((ok = gmx_stats_add_point(gstats, x[i], y[i],
187 (NULL((void*)0) != dx) ? dx[i] : 0,
188 (NULL((void*)0) != dy) ? dy[i] : 0)) != estatsOK)
189 {
190 return ok;
191 }
192 }
193 return estatsOK;
194}
195
196static int gmx_stats_compute(gmx_stats *stats, int weight)
197{
198 double yy, yx, xx, sx, sy, dy, chi2, chi2aa, d2;
199 double ssxx, ssyy, ssxy;
200 double w, wtot, yx_nw, sy_nw, sx_nw, yy_nw, xx_nw, dx2, dy2;
201 int i, N;
202
203 N = stats->np;
204 if (stats->computed == 0)
205 {
206 if (N < 1)
207 {
208 return estatsNO_POINTS;
209 }
210
211 xx = xx_nw = 0;
212 yy = yy_nw = 0;
213 yx = yx_nw = 0;
214 sx = sx_nw = 0;
215 sy = sy_nw = 0;
216 wtot = 0;
217 d2 = 0;
218 for (i = 0; (i < N); i++)
219 {
220 d2 += dsqr(stats->x[i]-stats->y[i]);
221 if ((stats->dy[i]) && (weight == elsqWEIGHT_Y))
222 {
223 w = 1/dsqr(stats->dy[i]);
224 }
225 else
226 {
227 w = 1;
228 }
229
230 wtot += w;
231
232 xx += w*dsqr(stats->x[i]);
233 xx_nw += dsqr(stats->x[i]);
234
235 yy += w*dsqr(stats->y[i]);
236 yy_nw += dsqr(stats->y[i]);
237
238 yx += w*stats->y[i]*stats->x[i];
239 yx_nw += stats->y[i]*stats->x[i];
240
241 sx += w*stats->x[i];
242 sx_nw += stats->x[i];
243
244 sy += w*stats->y[i];
245 sy_nw += stats->y[i];
246 }
247
248 /* Compute average, sigma and error */
249 stats->aver = sy_nw/N;
250 stats->sigma_aver = sqrt(yy_nw/N - dsqr(sy_nw/N));
251 stats->error = stats->sigma_aver/sqrt(N);
252
253 /* Compute RMSD between x and y */
254 stats->rmsd = sqrt(d2/N);
255
256 /* Correlation coefficient for data */
257 yx_nw /= N;
258 xx_nw /= N;
259 yy_nw /= N;
260 sx_nw /= N;
261 sy_nw /= N;
262 ssxx = N*(xx_nw - dsqr(sx_nw));
263 ssyy = N*(yy_nw - dsqr(sy_nw));
264 ssxy = N*(yx_nw - (sx_nw*sy_nw));
265 stats->Rdata = sqrt(dsqr(ssxy)/(ssxx*ssyy));
266
267 /* Compute straight line through datapoints, either with intercept
268 zero (result in aa) or with intercept variable (results in a
269 and b) */
270 yx = yx/wtot;
271 xx = xx/wtot;
272 sx = sx/wtot;
273 sy = sy/wtot;
274
275 stats->aa = (yx/xx);
276 stats->a = (yx-sx*sy)/(xx-sx*sx);
277 stats->b = (sy)-(stats->a)*(sx);
278
279 /* Compute chi2, deviation from a line y = ax+b. Also compute
280 chi2aa which returns the deviation from a line y = ax. */
281 chi2 = 0;
282 chi2aa = 0;
283 for (i = 0; (i < N); i++)
284 {
285 if (stats->dy[i] > 0)
286 {
287 dy = stats->dy[i];
288 }
289 else
290 {
291 dy = 1;
292 }
293 chi2aa += dsqr((stats->y[i]-(stats->aa*stats->x[i]))/dy);
294 chi2 += dsqr((stats->y[i]-(stats->a*stats->x[i]+stats->b))/dy);
295 }
296 if (N > 2)
297 {
298 stats->chi2 = sqrt(chi2/(N-2));
299 stats->chi2aa = sqrt(chi2aa/(N-2));
300
301 /* Look up equations! */
302 dx2 = (xx-sx*sx);
303 dy2 = (yy-sy*sy);
304 stats->sigma_a = sqrt(stats->chi2/((N-2)*dx2));
305 stats->sigma_b = stats->sigma_a*sqrt(xx);
306 stats->Rfit = fabs(ssxy)/sqrt(ssxx*ssyy);
307 /*stats->a*sqrt(dx2/dy2);*/
308 stats->Rfitaa = stats->aa*sqrt(dx2/dy2);
309 }
310 else
311 {
312 stats->chi2 = 0;
313 stats->chi2aa = 0;
314 stats->sigma_a = 0;
315 stats->sigma_b = 0;
316 stats->Rfit = 0;
317 stats->Rfitaa = 0;
318 }
319
320 stats->computed = 1;
321 }
322
323 return estatsOK;
324}
325
326int gmx_stats_get_ab(gmx_stats_t gstats, int weight,
327 real *a, real *b, real *da, real *db,
328 real *chi2, real *Rfit)
329{
330 gmx_stats *stats = (gmx_stats *) gstats;
331 int ok;
332
333 if ((ok = gmx_stats_compute(stats, weight)) != estatsOK)
334 {
335 return ok;
336 }
337 if (NULL((void*)0) != a)
338 {
339 *a = stats->a;
340 }
341 if (NULL((void*)0) != b)
342 {
343 *b = stats->b;
344 }
345 if (NULL((void*)0) != da)
346 {
347 *da = stats->sigma_a;
348 }
349 if (NULL((void*)0) != db)
350 {
351 *db = stats->sigma_b;
352 }
353 if (NULL((void*)0) != chi2)
354 {
355 *chi2 = stats->chi2;
356 }
357 if (NULL((void*)0) != Rfit)
358 {
359 *Rfit = stats->Rfit;
360 }
361
362 return estatsOK;
363}
364
365int gmx_stats_get_a(gmx_stats_t gstats, int weight, real *a, real *da,
366 real *chi2, real *Rfit)
367{
368 gmx_stats *stats = (gmx_stats *) gstats;
369 int ok;
370
371 if ((ok = gmx_stats_compute(stats, weight)) != estatsOK)
372 {
373 return ok;
374 }
375 if (NULL((void*)0) != a)
376 {
377 *a = stats->aa;
378 }
379 if (NULL((void*)0) != da)
380 {
381 *da = stats->sigma_aa;
382 }
383 if (NULL((void*)0) != chi2)
384 {
385 *chi2 = stats->chi2aa;
386 }
387 if (NULL((void*)0) != Rfit)
388 {
389 *Rfit = stats->Rfitaa;
390 }
391
392 return estatsOK;
393}
394
395int gmx_stats_get_average(gmx_stats_t gstats, real *aver)
396{
397 gmx_stats *stats = (gmx_stats *) gstats;
398 int ok;
399
400 if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
401 {
402 return ok;
403 }
404
405 *aver = stats->aver;
406
407 return estatsOK;
408}
409
410int gmx_stats_get_ase(gmx_stats_t gstats, real *aver, real *sigma, real *error)
411{
412 gmx_stats *stats = (gmx_stats *) gstats;
413 int ok;
414
415 if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
416 {
417 return ok;
418 }
419
420 if (NULL((void*)0) != aver)
421 {
422 *aver = stats->aver;
423 }
424 if (NULL((void*)0) != sigma)
425 {
426 *sigma = stats->sigma_aver;
427 }
428 if (NULL((void*)0) != error)
429 {
430 *error = stats->error;
431 }
432
433 return estatsOK;
434}
435
436int gmx_stats_get_sigma(gmx_stats_t gstats, real *sigma)
437{
438 gmx_stats *stats = (gmx_stats *) gstats;
439 int ok;
440
441 if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
442 {
443 return ok;
444 }
445
446 *sigma = stats->sigma_aver;
447
448 return estatsOK;
449}
450
451int gmx_stats_get_error(gmx_stats_t gstats, real *error)
452{
453 gmx_stats *stats = (gmx_stats *) gstats;
454 int ok;
455
456 if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
457 {
458 return ok;
459 }
460
461 *error = stats->error;
462
463 return estatsOK;
464}
465
466int gmx_stats_get_corr_coeff(gmx_stats_t gstats, real *R)
467{
468 gmx_stats *stats = (gmx_stats *) gstats;
469 int ok;
470
471 if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
472 {
473 return ok;
474 }
475
476 *R = stats->Rdata;
477
478 return estatsOK;
479}
480
481int gmx_stats_get_rmsd(gmx_stats_t gstats, real *rmsd)
482{
483 gmx_stats *stats = (gmx_stats *) gstats;
484 int ok;
485
486 if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
487 {
488 return ok;
489 }
490
491 *rmsd = stats->rmsd;
492
493 return estatsOK;
494}
495
496int gmx_stats_dump_xy(gmx_stats_t gstats, FILE *fp)
497{
498 gmx_stats *stats = (gmx_stats *) gstats;
499 int i, ok;
500
501 for (i = 0; (i < stats->np); i++)
502 {
503 fprintf(fp, "%12g %12g %12g %12g\n", stats->x[i], stats->y[i],
504 stats->dx[i], stats->dy[i]);
505 }
506
507 return estatsOK;
508}
509
510int gmx_stats_remove_outliers(gmx_stats_t gstats, double level)
511{
512 gmx_stats *stats = (gmx_stats *) gstats;
513 int i, iter = 1, done = 0, ok;
514 real rmsd, r;
515
516 while ((stats->np >= 10) && !done)
517 {
518 if ((ok = gmx_stats_get_rmsd(gstats, &rmsd)) != estatsOK)
519 {
520 return ok;
521 }
522 done = 1;
523 for (i = 0; (i < stats->np); )
524 {
525 r = fabs(stats->x[i]-stats->y[i]);
526 if (r > level*rmsd)
527 {
528 fprintf(stderrstderr, "Removing outlier, iter = %d, rmsd = %g, x = %g, y = %g\n",
529 iter, rmsd, stats->x[i], stats->y[i]);
530 if (i < stats->np-1)
531 {
532 stats->x[i] = stats->x[stats->np-1];
533 stats->y[i] = stats->y[stats->np-1];
534 stats->dx[i] = stats->dx[stats->np-1];
535 stats->dy[i] = stats->dy[stats->np-1];
536 }
537 stats->np--;
538 done = 0;
539 }
540 else
541 {
542 i++;
543 }
544 }
545 iter++;
546 }
547
548 return estatsOK;
549}
550
551int gmx_stats_make_histogram(gmx_stats_t gstats, real binwidth, int *nb,
552 int ehisto, int normalized, real **x, real **y)
553{
554 gmx_stats *stats = (gmx_stats *) gstats;
555 int i, ok, index = 0, nbins = *nb, *nindex;
556 double minx, maxx, maxy, miny, delta, dd, minh;
557
558 if (((binwidth <= 0) && (nbins <= 0)) ||
559 ((binwidth > 0) && (nbins > 0)))
560 {
561 return estatsINVALID_INPUT;
562 }
563 if (stats->np <= 2)
564 {
565 return estatsNO_POINTS;
566 }
567 minx = maxx = stats->x[0];
568 miny = maxy = stats->y[0];
569 for (i = 1; (i < stats->np); i++)
570 {
571 miny = (stats->y[i] < miny) ? stats->y[i] : miny;
572 maxy = (stats->y[i] > maxy) ? stats->y[i] : maxy;
573 minx = (stats->x[i] < minx) ? stats->x[i] : minx;
574 maxx = (stats->x[i] > maxx) ? stats->x[i] : maxx;
575 }
576 if (ehisto == ehistoX)
577 {
578 delta = maxx-minx;
579 minh = minx;
580 }
581 else if (ehisto == ehistoY)
582 {
583 delta = maxy-miny;
584 minh = miny;
585 }
586 else
587 {
588 return estatsINVALID_INPUT;
589 }
590
591 if (binwidth == 0)
592 {
593 binwidth = (delta)/nbins;
594 }
595 else
596 {
597 nbins = gmx_dnint((delta)/binwidth + 0.5);
598 }
599 snew(*x, nbins)(*x) = save_calloc("*x", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 599, (nbins), sizeof(*(*x)));
600 snew(nindex, nbins)(nindex) = save_calloc("nindex", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 600, (nbins), sizeof(*(nindex)));
601 for (i = 0; (i < nbins); i++)
602 {
603 (*x)[i] = minh + binwidth*(i+0.5);
604 }
605 if (normalized == 0)
606 {
607 dd = 1;
608 }
609 else
610 {
611 dd = 1.0/(binwidth*stats->np);
612 }
613
614 snew(*y, nbins)(*y) = save_calloc("*y", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 614, (nbins), sizeof(*(*y)));
615 for (i = 0; (i < stats->np); i++)
616 {
617 if (ehisto == ehistoY)
618 {
619 index = (stats->y[i]-miny)/binwidth;
620 }
621 else if (ehisto == ehistoX)
622 {
623 index = (stats->x[i]-minx)/binwidth;
624 }
625 if (index < 0)
626 {
627 index = 0;
628 }
629 if (index > nbins-1)
630 {
631 index = nbins-1;
632 }
633 (*y)[index] += dd;
634 nindex[index]++;
635 }
636 if (*nb == 0)
637 {
638 *nb = nbins;
639 }
640 for (i = 0; (i < nbins); i++)
641 {
642 if (nindex[i] > 0)
643 {
644 (*y)[i] /= nindex[i];
645 }
646 }
647
648 sfree(nindex)save_free("nindex", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 648, (nindex));
649
650 return estatsOK;
651}
652
653static const char *stats_error[estatsNR] =
654{
655 "All well in STATS land",
656 "No points",
657 "Not enough memory",
658 "Invalid histogram input",
659 "Unknown error",
660 "Not implemented yet"
661};
662
663const char *gmx_stats_message(int estats)
664{
665 if ((estats >= 0) && (estats < estatsNR))
666 {
667 return stats_error[estats];
668 }
669 else
670 {
671 return stats_error[estatsERROR];
672 }
673}
674
675/* Old convenience functions, should be merged with the core
676 statistics above. */
677int lsq_y_ax(int n, real x[], real y[], real *a)
678{
679 gmx_stats_t lsq = gmx_stats_init();
680 int ok;
681 real da, chi2, Rfit;
682
683 gmx_stats_add_points(lsq, n, x, y, 0, 0);
684 if ((ok = gmx_stats_get_a(lsq, elsqWEIGHT_NONE, a, &da, &chi2, &Rfit)) != estatsOK)
685 {
686 return ok;
687 }
688
689 /* int i;
690 double xx,yx;
691
692 yx=xx=0.0;
693 for (i=0; i<n; i++) {
694 yx+=y[i]*x[i];
695 xx+=x[i]*x[i];
696 }
697 * a=yx/xx;
698 */
699 return estatsOK;
700}
701
702static int low_lsq_y_ax_b(int n, real *xr, double *xd, real yr[],
703 real *a, real *b, real *r, real *chi2)
704{
705 int i, ok;
706 gmx_stats_t lsq;
707
708 lsq = gmx_stats_init();
709 for (i = 0; (i < n); i++)
3
←
Assuming 'i' is < 'n'
→
4
←
Loop condition is true. Entering loop body
→
710 {
711 if ((ok = gmx_stats_add_point(lsq, (NULL((void*)0) != xd) ? xd[i] : xr[i], yr[i], 0, 0))
5
←
'?' condition is false
→
6
←
Array access (from variable 'xr') results in a null pointer dereference
712 != estatsOK)
713 {
714 return ok;
715 }
716 }
717 if ((ok = gmx_stats_get_ab(lsq, elsqWEIGHT_NONE, a, b, NULL((void*)0), NULL((void*)0), chi2, r)) != estatsOK)
718 {
719 return ok;
720 }
721
722 return estatsOK;
723 /*
724 double x,y,yx,xx,yy,sx,sy,chi2;
725
726 yx=xx=yy=sx=sy=0.0;
727 for (i=0; i<n; i++) {
728 if (xd != NULL) {
729 x = xd[i];
730 } else {
731 x = xr[i];
732 }
733 y = yr[i];
734
735 yx += y*x;
736 xx += x*x;
737 yy += y*y;
738 sx += x;
739 sy += y;
740 }
741 * a = (n*yx-sy*sx)/(n*xx-sx*sx);
742 * b = (sy-(*a)*sx)/n;
743 * r = sqrt((xx-sx*sx)/(yy-sy*sy));
744
745 chi2 = 0;
746 if (xd != NULL) {
747 for(i=0; i<n; i++)
748 chi2 += dsqr(yr[i] - ((*a)*xd[i] + (*b)));
749 } else {
750 for(i=0; i<n; i++)
751 chi2 += dsqr(yr[i] - ((*a)*xr[i] + (*b)));
752 }
753
754 if (n > 2)
755 return sqrt(chi2/(n-2));
756 else
757 return 0;
758 */
759}
760
761int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b, real *r, real *chi2)
762{
763 return low_lsq_y_ax_b(n, x, NULL((void*)0), y, a, b, r, chi2);
764}
765
766int lsq_y_ax_b_xdouble(int n, double x[], real y[], real *a, real *b,
767 real *r, real *chi2)
768{
769 return low_lsq_y_ax_b(n, NULL((void*)0), x, y, a, b, r, chi2);
1
Within the expansion of the macro 'NULL':
→
a
Passing null pointer value via 2nd parameter 'xr'
2
←
Calling 'low_lsq_y_ax_b'
→
770}
771
772int lsq_y_ax_b_error(int n, real x[], real y[], real dy[],
773 real *a, real *b, real *da, real *db,
774 real *r, real *chi2)
775{
776 gmx_stats_t lsq;
777 int i, ok;
778
779 lsq = gmx_stats_init();
780 for (i = 0; (i < n); i++)
781 {
782 if ((ok = gmx_stats_add_point(lsq, x[i], y[i], 0, dy[i])) != estatsOK)
783 {
784 return ok;
785 }
786 }
787 if ((ok = gmx_stats_get_ab(lsq, elsqWEIGHT_Y, a, b, da, db, chi2, r)) != estatsOK)
788 {
789 return ok;
790 }
791 if ((ok = gmx_stats_done(lsq)) != estatsOK)
792 {
793 return ok;
794 }
795 sfree(lsq)save_free("lsq", "/home/alexxy/Develop/gromacs/src/gromacs/statistics/statistics.c"
, 795, (lsq));
796
797 return estatsOK;
798 /*
799 double sxy,sxx,syy,sx,sy,w,s_2,dx2,dy2,mins;
800
801 sxy=sxx=syy=sx=sy=w=0.0;
802 mins = dy[0];
803 for(i=1; (i<n); i++)
804 mins = min(mins,dy[i]);
805 if (mins <= 0)
806 gmx_fatal(FARGS,"Zero or negative weigths in linear regression analysis");
807
808 for (i=0; i<n; i++) {
809 s_2 = dsqr(1.0/dy[i]);
810 sxx += s_2*dsqr(x[i]);
811 sxy += s_2*y[i]*x[i];
812 syy += s_2*dsqr(y[i]);
813 sx += s_2*x[i];
814 sy += s_2*y[i];
815 w += s_2;
816 }
817 sxx = sxx/w;
818 sxy = sxy/w;
819 syy = syy/w;
820 sx = sx/w;
821 sy = sy/w;
822 dx2 = (sxx-sx*sx);
823 dy2 = (syy-sy*sy);
824 * a=(sxy-sy*sx)/dx2;
825 * b=(sy-(*a)*sx);
826
827 * chi2=0;
828 for(i=0; i<n; i++)
829 * chi2+=dsqr((y[i]-((*a)*x[i]+(*b)))/dy[i]);
830 * chi2 = *chi2/w;
831
832 * da = sqrt(*chi2/((n-2)*dx2));
833 * db = *da*sqrt(sxx);
834 * r = *a*sqrt(dx2/dy2);
835
836 if (debug)
837 fprintf(debug,"sx = %g, sy = %g, sxy = %g, sxx = %g, w = %g\n"
838 "chi2 = %g, dx2 = %g\n",
839 sx,sy,sxy,sxx,w,*chi2,dx2);
840
841 if (n > 2)
842 * chi2 = sqrt(*chi2/(n-2));
843 else
844 * chi2 = 0;
845 */
846}