1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 This file contains general purpose routines for input, manipulation, and
10 output of probability density functions.
11
12 File: pdf.c
13 Author: B. Douglas Ward
14 Date: 28 January 2000
15 */
16
17 #ifndef USE_QUIET
18 # define quiet 0
19 #endif
20
21
22 /*---------------------------------------------------------------------------*/
23 /*
24 Routine to print an error message and stop.
25 */
26
PDF_error(char * message)27 void PDF_error (char * message)
28 {
29 printf ("PDF error: %s \n", message);
30 exit (1);
31 }
32
33
34 /*---------------------------------------------------------------------------*/
35 /*
36 Initialize pdf data structure.
37 */
38
PDF_initialize(pdf * p)39 void PDF_initialize (pdf * p)
40 {
41 p->nbin = 0;
42 p->prob = NULL;
43 p->lower_bnd = 0.0;
44 p->upper_bnd = 0.0;
45 p->width = 0.0;
46
47 return;
48 }
49
50
51 /*---------------------------------------------------------------------------*/
52 /*
53 Destroy pdf data structure by deallocating memory.
54 */
55
PDF_destroy(pdf * p)56 void PDF_destroy (pdf * p)
57 {
58
59 if (p->prob != NULL) free (p->prob);
60
61 PDF_initialize (p);
62
63 return;
64 }
65
66
67 /*---------------------------------------------------------------------------*/
68 /*
69 Normalize pdf to have unity area.
70 */
71
PDF_normalize(pdf * p)72 void PDF_normalize (pdf * p)
73 {
74 int ibin;
75 double sum;
76
77
78 sum = 0.0;
79
80 for (ibin = 0; ibin < p->nbin; ibin++)
81 sum += p->prob[ibin];
82
83
84 for (ibin = 0; ibin < p->nbin; ibin++)
85 p->prob[ibin] /= sum;
86
87 return;
88 }
89
90
91 /*---------------------------------------------------------------------------*/
92 /*
93 Create pdf data structure by allocating memory and initializing values.
94 */
95
PDF_create(int nbin,float * prob,float lower_bnd,float upper_bnd,pdf * p)96 void PDF_create (int nbin, float * prob, float lower_bnd, float upper_bnd,
97 pdf * p)
98 {
99 int ibin;
100
101
102 PDF_destroy (p);
103
104 p->nbin = nbin;
105
106 p->prob = (float *) malloc (sizeof(float) * nbin);
107 PDF_MTEST (p->prob);
108 for (ibin = 0; ibin < nbin; ibin++)
109 p->prob[ibin] = prob[ibin];
110
111 p->lower_bnd = lower_bnd;
112 p->upper_bnd = upper_bnd;
113
114 p->width = (upper_bnd - lower_bnd) / (nbin-1);
115
116 PDF_normalize (p);
117
118 return;
119 }
120
121
122 /*---------------------------------------------------------------------------*/
123 /*
124 Copy pdf data structure from p to pc.
125 */
126
PDF_copy(pdf p,pdf * pc)127 void PDF_copy (pdf p, pdf * pc)
128 {
129 PDF_create (p.nbin, p.prob, p.lower_bnd, p.upper_bnd, pc);
130
131 return;
132 }
133
134
135 /*---------------------------------------------------------------------------*/
136 /*
137 Convert bin number to value of independent variable at center of the bin.
138 */
139
PDF_ibin_to_xvalue(pdf p,int ibin)140 float PDF_ibin_to_xvalue (pdf p, int ibin)
141 {
142 float xvalue;
143
144 xvalue = p.lower_bnd + ibin*p.width;
145
146 return (xvalue);
147 }
148
149
150 /*---------------------------------------------------------------------------*/
151 /*
152 Convert value of independent variable to bin number.
153 */
154
PDF_xvalue_to_ibin(pdf p,float xvalue)155 int PDF_xvalue_to_ibin (pdf p, float xvalue)
156 {
157 int ibin;
158
159 ibin = floor ( (xvalue - p.lower_bnd) / p.width + 0.5);
160
161 return (ibin);
162 }
163
164
165 /*---------------------------------------------------------------------------*/
166 /*
167 Convert value of independent variable to probability.
168 */
169
PDF_xvalue_to_pvalue(pdf p,float xvalue)170 float PDF_xvalue_to_pvalue (pdf p, float xvalue)
171 {
172 int ibin;
173 float pvalue;
174
175 ibin = PDF_xvalue_to_ibin (p, xvalue);
176
177 if ((ibin < 0) || (ibin >= p.nbin))
178 pvalue = 0.0;
179 else
180 pvalue = p.prob[ibin];
181
182 return (pvalue);
183 }
184
185
186 /*---------------------------------------------------------------------------*/
187 /*
188 Print contents of pdf p to screen.
189 */
190
PDF_print(pdf p)191 void PDF_print (pdf p)
192 {
193 int ibin;
194
195
196 if( !quiet ){
197 printf ("Number of bins = %d \n", p.nbin);
198 printf ("Lower bound = %f \n", p.lower_bnd);
199 printf ("Upper bound = %f \n", p.upper_bnd);
200 printf ("Bin width = %f \n", p.width);
201
202 /*
203 printf ("%3s %10.6s %10.6s \n", "i", "x[i]", "p[i]");
204
205 for (ibin = 0; ibin < p.nbin; ibin++)
206 printf ("%3d %10.6f %10.6f \n",
207 ibin, PDF_ibin_to_xvalue(p, ibin), p.prob[ibin]);
208 */
209 }
210
211 return;
212 }
213
214
215 /*---------------------------------------------------------------------------*/
216 /*
217 Print contents of string str and pdf p to screen.
218 */
219
PDF_sprint(char * str,pdf p)220 void PDF_sprint (char * str, pdf p)
221 {
222 if( quiet ) return ;
223 printf ("%s \n", str);
224
225 PDF_print (p);
226
227 return;
228 }
229
230
231 /*---------------------------------------------------------------------------*/
232 /*
233 Write contents of pdf p to specified file.
234 */
235
PDF_write_file(char * filename,pdf p)236 void PDF_write_file (char * filename, pdf p)
237 {
238 int ibin;
239 FILE * outfile = NULL;
240
241
242 outfile = fopen (filename, "w");
243 if (!outfile) {
244 fprintf (stderr, "\n"
245 "*****************************\n"
246 "Error:\n"
247 "Failed to open %s for output.\n"
248 "Check for write permissions.\n"
249 "*****************************\n"
250 "\n", filename);
251 return;
252 }
253 for (ibin = 0; ibin < p.nbin; ibin++)
254 fprintf (outfile, "%d %f %f \n",
255 ibin, PDF_ibin_to_xvalue(p, ibin), p.prob[ibin]);
256
257
258 fclose (outfile);
259
260 return;
261 }
262
263
264 /*---------------------------------------------------------------------------*/
265 /*
266 Simple smoothing of the pdf estimate.
267 */
268
PDF_smooth(pdf * p)269 void PDF_smooth (pdf * p)
270 {
271 float * sprob;
272 int ibin;
273
274
275 sprob = (float *) malloc (sizeof(float) * p->nbin);
276
277 sprob[0] = 0.5*(p->prob[0] + p->prob[1]);
278 sprob[p->nbin-1] = 0.5*(p->prob[p->nbin-2] + p->prob[p->nbin-1]);
279
280 for (ibin = 1; ibin < p->nbin-1; ibin++)
281 sprob[ibin] = 0.25 * (p->prob[ibin-1] + 2*p->prob[ibin] + p->prob[ibin+1]);
282
283 free (p->prob);
284 p->prob = sprob;
285
286 PDF_normalize (p);
287
288 return;
289 }
290
291
292 /*---------------------------------------------------------------------------*/
293 /*
294 Trim the pdf by removing the extreme lower and extreme upper values.
295 */
296
PDF_trim(float lower_per,float upper_per,pdf * p)297 void PDF_trim (float lower_per, float upper_per, pdf * p)
298 {
299 int ibin;
300 float * fbin = NULL;
301 float cum_prob;
302 float lower_bnd, upper_bnd;
303 int lo_bin=0, hi_bin=0;
304
305
306 /*----- Trim lower values -----*/
307 cum_prob = 0.0;
308 for (ibin = 0; ibin < p->nbin; ibin++)
309 {
310 cum_prob += p->prob[ibin];
311 p->prob[ibin] = 0.0;
312 if (cum_prob > lower_per)
313 {
314 lo_bin = ibin + 1;
315 break;
316 }
317 }
318
319
320 /*----- Trim upper values -----*/
321 cum_prob = 0.0;
322 for (ibin = p->nbin-1; ibin >= 0; ibin--)
323 {
324 cum_prob += p->prob[ibin];
325 p->prob[ibin] = 0.0;
326 if (cum_prob > 1.0 - upper_per)
327 {
328 hi_bin = ibin - 1;
329 break;
330 }
331 }
332
333
334 /*----- Reset lower and upper bounds -----*/
335 lower_bnd = PDF_ibin_to_xvalue (*p, lo_bin);
336 upper_bnd = PDF_ibin_to_xvalue (*p, hi_bin);
337
338 p->lower_bnd = lower_bnd;
339 p->upper_bnd = upper_bnd;
340 p->nbin = hi_bin - lo_bin + 1;
341
342
343 /*----- Copy data -----*/
344 fbin = (float *) malloc (sizeof(float) * p->nbin);
345 for (ibin = 0; ibin < p->nbin; ibin++)
346 fbin[ibin] = p->prob[ibin+lo_bin];
347
348
349 /*----- Reset pointer to data -----*/
350 free (p->prob);
351 p->prob = fbin;
352
353
354 /*----- Normalize to unity area -----*/
355 PDF_normalize (p);
356
357
358 return;
359 }
360
361
362 /*---------------------------------------------------------------------------*/
363 /*
364 Determine the range of values in the input short array.
365 */
366
PDF_short_range(int npts,short * sarray,short * min_val,short * max_val)367 void PDF_short_range (int npts, short * sarray,
368 short * min_val, short * max_val)
369 {
370 int ipt;
371
372
373 *min_val = sarray[0];
374 *max_val = sarray[0];
375
376 for (ipt = 1; ipt < npts; ipt++)
377 {
378 if (sarray[ipt] < *min_val) *min_val = sarray[ipt];
379 if (sarray[ipt] > *max_val) *max_val = sarray[ipt];
380 }
381
382 return;
383 }
384
385
386 /*---------------------------------------------------------------------------*/
387 /*
388 Determine the range of values in the input float array.
389 */
390
PDF_float_range(int npts,float * farray,float * min_val,float * max_val)391 void PDF_float_range (int npts, float * farray,
392 float * min_val, float * max_val)
393 {
394 int ipt;
395
396
397 *min_val = farray[0];
398 *max_val = farray[0];
399
400 for (ipt = 1; ipt < npts; ipt++)
401 {
402 if (farray[ipt] < *min_val) *min_val = farray[ipt];
403 if (farray[ipt] > *max_val) *max_val = farray[ipt];
404 }
405
406 return;
407 }
408
409
410 /*---------------------------------------------------------------------------*/
411 /*
412 Estimate the pdf corresponding to the input short array.
413 */
414
PDF_short_to_pdf(int npts,short * sarray,pdf * p)415 void PDF_short_to_pdf (int npts, short * sarray, pdf * p)
416 {
417 const int MIN_COUNT = 5;
418 const int MIN_BINS = 5;
419 int ipt, ibin, count;
420 float * fbin = NULL;
421 int num_bins;
422 short lower_lim, upper_lim;
423 char message[80];
424
425
426 /*----- Make histogram of input short array -----*/
427 PDF_short_range (npts, sarray, &lower_lim, &upper_lim);
428 num_bins = upper_lim - lower_lim + 1;
429 if (num_bins < MIN_BINS)
430 {
431 sprintf (message, "histogram contains only %d bins", num_bins);
432 PDF_error (message);
433 }
434
435 fbin = (float *) malloc (sizeof(float) * num_bins); PDF_MTEST (fbin);
436 for (ibin = 0; ibin < num_bins; ibin++)
437 fbin[ibin] = 0.0;
438
439 count = 0;
440 for (ipt = 0; ipt < npts; ipt++)
441 {
442 ibin = sarray[ipt] - lower_lim;
443 if ((ibin >= 0) && (ibin < num_bins))
444 {
445 fbin[ibin] += 1.0;
446 count++;
447 }
448 }
449
450
451 /*----- Check for too few points -----*/
452 if (count < MIN_COUNT)
453 {
454 sprintf (message, "histogram contains only %d points", count);
455 PDF_error (message);
456 }
457
458
459 /*----- Create PDF -----*/
460 PDF_create (num_bins, fbin, (float) lower_lim, (float) upper_lim, p);
461
462
463 /*----- Release memory -----*/
464 free (fbin); fbin = NULL;
465
466
467 return;
468 }
469
470
471 /*---------------------------------------------------------------------------*/
472 /*
473 Estimate the pdf corresponding to the input float array.
474 */
PDF_float_to_pdf(int npts,float * farray,int num_bins,pdf * p)475 void PDF_float_to_pdf (int npts, float * farray, int num_bins, pdf * p)
476 {
477 const int MIN_COUNT = 5;
478 const int MIN_BINS = 5;
479 int ipt, ibin, count;
480 float * fbin = NULL;
481 float width;
482 float min_val, max_val;
483 char message[80];
484
485
486 /*----- Make histogram of input float array -----*/
487 if (num_bins < MIN_BINS)
488 {
489 sprintf (message, "histogram contains only %d bins", num_bins);
490 PDF_error (message);
491 }
492
493 fbin = (float *) malloc (sizeof(float) * num_bins); PDF_MTEST (fbin);
494 for (ibin = 0; ibin < num_bins; ibin++)
495 fbin[ibin] = 0.0;
496
497 PDF_float_range (npts, farray, &min_val, &max_val);
498 width = (max_val - min_val) / num_bins;
499
500 count = 0;
501 for (ipt = 0; ipt < npts; ipt++)
502 {
503 ibin = (farray[ipt] - min_val) / width;
504 if ((ibin >= 0) && (ibin < num_bins))
505 {
506 fbin[ibin] += 1.0;
507 count++;
508 }
509 }
510
511
512 /*----- Check for too few points -----*/
513 if (count < MIN_COUNT)
514 {
515 sprintf (message, "histogram contains only %d points", count);
516 PDF_error (message);
517 }
518
519
520 /*----- Create PDF -----*/
521 PDF_create (num_bins, fbin, min_val, max_val, p);
522
523 /*----- Release memory -----*/
524 free (fbin); fbin = NULL;
525
526
527 return;
528 }
529
530
531 /*---------------------------------------------------------------------------*/
532 /*
533 Find extrema of pdf function.
534 */
535
PDF_find_extrema(pdf p,int * num_min,int * pdf_min,int * num_max,int * pdf_max)536 void PDF_find_extrema (pdf p, int * num_min, int * pdf_min,
537 int * num_max, int * pdf_max)
538 {
539 int ibin;
540 int i;
541
542
543 *num_min = 0;
544 *num_max = 0;
545
546 for (ibin = 1; ibin < p.nbin-1; ibin++)
547 {
548 if ((p.prob[ibin] < p.prob[ibin-1]) && (p.prob[ibin] < p.prob[ibin+1]))
549 {
550 pdf_min[*num_min] = ibin;
551 (*num_min)++;
552 }
553
554 if ((p.prob[ibin] > p.prob[ibin-1]) && (p.prob[ibin] > p.prob[ibin+1]))
555 {
556 pdf_max[*num_max] = ibin;
557 (*num_max)++;
558 }
559 }
560
561 if( !quiet ){
562 printf ("\nExtrema of PDF: \n");
563 printf ("\nNum Local Min = %d \n", *num_min);
564 for (i = 0; i < *num_min; i++)
565 {
566 ibin = pdf_min[i];
567 printf ("x[%3d] = %8.3f p[%3d] = %12.6f \n",
568 ibin, PDF_ibin_to_xvalue(p, ibin), ibin, p.prob[ibin]);
569 }
570
571 printf ("\nNum Local Max = %d \n", *num_max);
572 for (i = 0; i < *num_max; i++)
573 {
574 ibin = pdf_max[i];
575 printf ("x[%3d] = %8.3f p[%3d] = %12.6f \n",
576 ibin, PDF_ibin_to_xvalue(p, ibin), ibin, p.prob[ibin]);
577 }
578 }
579
580 }
581
582
583 /*---------------------------------------------------------------------------*/
584 /*
585 Find bimodality of pdf function (if possible).
586 */
587
PDF_find_bimodal(pdf p,int * gmax,int * wmax)588 int PDF_find_bimodal (pdf p, int * gmax, int * wmax)
589 {
590 const int NPTS = 12;
591 int * pdf_min = NULL, * pdf_max = NULL;
592 int num_min, num_max;
593 int imax, temp;
594
595
596 pdf_min = (int *) malloc (sizeof(int) * p.nbin);
597 pdf_max = (int *) malloc (sizeof(int) * p.nbin);
598
599 PDF_find_extrema (p, &num_min, pdf_min, &num_max, pdf_max);
600
601
602 if (num_max >= 2)
603 {
604 if (p.prob[pdf_max[1]] >= p.prob[pdf_max[0]])
605 {
606 *wmax = pdf_max[1];
607 *gmax = pdf_max[0];
608 }
609 else
610 {
611 *wmax = pdf_max[0];
612 *gmax = pdf_max[1];
613 }
614
615 if (num_max > 2)
616 {
617 for (imax = 2; imax < num_max; imax++)
618 {
619 if (p.prob[pdf_max[imax]] >= p.prob[*wmax])
620 {
621 *gmax = *wmax;
622 *wmax = pdf_max[imax];
623 }
624 else if (p.prob[pdf_max[imax]] >= p.prob[*gmax])
625 {
626 *gmax = pdf_max[imax];
627 }
628 }
629 }
630
631 if (*gmax > *wmax)
632 {
633 temp = *gmax;
634 *gmax = *wmax;
635 *wmax = temp;
636 }
637
638 } /* end if (num_max >= 2) */
639
640
641 free (pdf_min); pdf_min = NULL;
642 free (pdf_max); pdf_max = NULL;
643
644
645 if (num_max < 2) return (0);
646 else return (1);
647
648 }
649
650
651 /*---------------------------------------------------------------------------*/
652
653
654
655
656
657
658