00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This file contains general purpose routines for input, manipulation, and
00010   output of probability density functions.
00011 
00012   File:    pdf.c
00013   Author:  B. Douglas Ward
00014   Date:    28 January 2000
00015 */
00016 
00017 #ifndef USE_QUIET
00018 # define quiet 0
00019 #endif
00020 
00021 
00022 /*---------------------------------------------------------------------------*/
00023 /*
00024   Routine to print an error message and stop.
00025 */
00026 
00027 void PDF_error (char * message)
00028 {
00029   printf ("PDF error: %s \n", message);
00030   exit (1);
00031 }
00032 
00033 
00034 /*---------------------------------------------------------------------------*/
00035 /*
00036   Initialize pdf data structure.
00037 */
00038 
00039 void PDF_initialize (pdf * p)
00040 {
00041   p->nbin = 0;
00042   p->prob = NULL;
00043   p->lower_bnd = 0.0;
00044   p->upper_bnd = 0.0;
00045   p->width = 0.0;
00046 
00047   return;
00048 }
00049 
00050 
00051 /*---------------------------------------------------------------------------*/
00052 /*
00053   Destroy pdf data structure by deallocating memory.
00054 */
00055 
00056 void PDF_destroy (pdf * p)
00057 {
00058 
00059   if (p->prob != NULL)   free (p->prob);
00060 
00061   PDF_initialize (p);
00062 
00063   return;
00064 }
00065 
00066 
00067 /*---------------------------------------------------------------------------*/
00068 /*
00069   Normalize pdf to have unity area.
00070 */
00071 
00072 void PDF_normalize (pdf * p)
00073 {
00074   int ibin;
00075   double sum;
00076 
00077 
00078   sum = 0.0;
00079 
00080   for (ibin = 0;  ibin < p->nbin;  ibin++)
00081     sum += p->prob[ibin];
00082 
00083 
00084   for (ibin = 0;  ibin < p->nbin;  ibin++)
00085     p->prob[ibin] /= sum;
00086 
00087   return;
00088 }
00089 
00090 
00091 /*---------------------------------------------------------------------------*/
00092 /*
00093   Create pdf data structure by allocating memory and initializing values.
00094 */
00095 
00096 void PDF_create (int nbin, float * prob, float lower_bnd, float upper_bnd, 
00097                  pdf * p)
00098 {
00099   int ibin;
00100 
00101 
00102   PDF_destroy (p);
00103 
00104   p->nbin = nbin;
00105   
00106   p->prob = (float *) malloc (sizeof(float) * nbin);
00107   PDF_MTEST (p->prob);
00108   for (ibin = 0;  ibin < nbin;  ibin++)
00109     p->prob[ibin] = prob[ibin];
00110 
00111   p->lower_bnd = lower_bnd;
00112   p->upper_bnd = upper_bnd;
00113 
00114   p->width = (upper_bnd - lower_bnd) / (nbin-1);
00115 
00116   PDF_normalize (p);
00117   
00118   return;
00119 }
00120 
00121 
00122 /*---------------------------------------------------------------------------*/
00123 /*
00124   Copy pdf data structure from p to pc.
00125 */
00126 
00127 void PDF_copy (pdf p, pdf * pc)
00128 {
00129   PDF_create (p.nbin, p.prob, p.lower_bnd, p.upper_bnd, pc);
00130 
00131   return;
00132 }
00133 
00134 
00135 /*---------------------------------------------------------------------------*/
00136 /*
00137   Convert bin number to value of independent variable at center of the bin.
00138 */
00139 
00140 float PDF_ibin_to_xvalue (pdf p, int ibin)
00141 {
00142   float xvalue;
00143 
00144   xvalue = p.lower_bnd + ibin*p.width;
00145 
00146   return (xvalue);
00147 }
00148 
00149   
00150 /*---------------------------------------------------------------------------*/
00151 /*
00152   Convert value of independent variable to bin number.
00153 */
00154 
00155 int PDF_xvalue_to_ibin (pdf p, float xvalue)
00156 {
00157   int ibin;
00158 
00159   ibin = floor ( (xvalue - p.lower_bnd) / p.width + 0.5);
00160 
00161   return (ibin);
00162 }
00163 
00164   
00165 /*---------------------------------------------------------------------------*/
00166 /*
00167   Convert value of independent variable to probability.
00168 */
00169 
00170 float PDF_xvalue_to_pvalue (pdf p, float xvalue)
00171 {
00172   int ibin;
00173   float pvalue;
00174 
00175   ibin = PDF_xvalue_to_ibin (p, xvalue);
00176 
00177   if ((ibin < 0) || (ibin >= p.nbin))
00178     pvalue = 0.0;
00179   else
00180     pvalue = p.prob[ibin];
00181 
00182   return (pvalue);
00183 }
00184 
00185   
00186 /*---------------------------------------------------------------------------*/
00187 /*
00188   Print contents of pdf p to screen.
00189 */
00190 
00191 void PDF_print (pdf p)
00192 {
00193   int ibin;
00194 
00195 
00196   if( !quiet ){
00197    printf ("Number of bins = %d \n", p.nbin);
00198    printf ("Lower bound    = %f \n", p.lower_bnd);
00199    printf ("Upper bound    = %f \n", p.upper_bnd);
00200    printf ("Bin width      = %f \n", p.width);
00201   
00202    /*
00203    printf ("%3s   %10.6s   %10.6s \n", "i", "x[i]", "p[i]");
00204  
00205    for (ibin = 0;  ibin < p.nbin;  ibin++)
00206      printf ("%3d   %10.6f   %10.6f \n", 
00207              ibin, PDF_ibin_to_xvalue(p, ibin), p.prob[ibin]);
00208    */
00209   }
00210 
00211   return;
00212 }
00213 
00214 
00215 /*---------------------------------------------------------------------------*/
00216 /*
00217   Print contents of string str and pdf p to screen.
00218 */
00219 
00220 void PDF_sprint (char * str, pdf p)
00221 {
00222   if( quiet ) return ;
00223   printf ("%s \n", str);
00224 
00225   PDF_print (p);  
00226 
00227   return;
00228 }
00229 
00230 
00231 /*---------------------------------------------------------------------------*/
00232 /*
00233   Write contents of pdf p to specified file.
00234 */
00235 
00236 void PDF_write_file (char * filename, pdf p)
00237 {
00238   int ibin;
00239   FILE * outfile = NULL;
00240 
00241 
00242   outfile = fopen (filename, "w");
00243 
00244   for (ibin = 0;  ibin < p.nbin;  ibin++)
00245     fprintf (outfile, "%d  %f  %f \n", 
00246             ibin, PDF_ibin_to_xvalue(p, ibin), p.prob[ibin]);
00247   
00248 
00249   fclose (outfile);
00250 
00251   return;
00252 }
00253 
00254 
00255 /*---------------------------------------------------------------------------*/
00256 /*
00257   Simple smoothing of the pdf estimate.
00258 */
00259 
00260 void PDF_smooth (pdf * p)
00261 {
00262   float * sprob;
00263   int ibin;
00264 
00265 
00266   sprob = (float *) malloc (sizeof(float) * p->nbin);
00267 
00268   sprob[0] = 0.5*(p->prob[0] + p->prob[1]);
00269   sprob[p->nbin-1] = 0.5*(p->prob[p->nbin-2] + p->prob[p->nbin-1]);
00270 
00271   for (ibin = 1;  ibin < p->nbin-1;  ibin++)
00272     sprob[ibin] = 0.25 * (p->prob[ibin-1] + 2*p->prob[ibin] + p->prob[ibin+1]);
00273   
00274   free (p->prob);
00275   p->prob = sprob;
00276 
00277   PDF_normalize (p);
00278 
00279   return;
00280 }
00281 
00282 
00283 /*---------------------------------------------------------------------------*/
00284 /*
00285   Trim the pdf by removing the extreme lower and extreme upper values.
00286 */
00287 
00288 void PDF_trim (float lower_per, float upper_per, pdf * p)
00289 {
00290   int ibin;
00291   float * fbin = NULL;
00292   float cum_prob;
00293   float lower_bnd, upper_bnd;
00294   int lo_bin, hi_bin;
00295 
00296 
00297   /*----- Trim lower values -----*/
00298   cum_prob = 0.0;
00299   for (ibin = 0;  ibin < p->nbin;  ibin++)
00300     {
00301       cum_prob += p->prob[ibin];
00302       p->prob[ibin] = 0.0;
00303       if (cum_prob > lower_per)
00304         {
00305           lo_bin = ibin + 1;
00306           break;
00307         }
00308     }
00309 
00310 
00311   /*----- Trim upper values -----*/
00312   cum_prob = 0.0;
00313   for (ibin = p->nbin-1;  ibin >= 0;  ibin--)
00314     {
00315       cum_prob += p->prob[ibin];
00316       p->prob[ibin] = 0.0;
00317       if (cum_prob > 1.0 - upper_per)
00318         {
00319           hi_bin = ibin - 1;
00320           break;
00321         }
00322     }
00323 
00324   
00325   /*----- Reset lower and upper bounds -----*/  
00326   lower_bnd = PDF_ibin_to_xvalue (*p, lo_bin);
00327   upper_bnd = PDF_ibin_to_xvalue (*p, hi_bin);
00328 
00329   p->lower_bnd = lower_bnd;
00330   p->upper_bnd = upper_bnd;
00331   p->nbin = hi_bin - lo_bin + 1;
00332 
00333 
00334   /*----- Copy data -----*/
00335   fbin = (float *) malloc (sizeof(float) * p->nbin);
00336   for (ibin = 0;  ibin < p->nbin;  ibin++)
00337     fbin[ibin] = p->prob[ibin+lo_bin];
00338 
00339 
00340   /*----- Reset pointer to data -----*/
00341   free (p->prob);
00342   p->prob = fbin;
00343 
00344 
00345   /*----- Normalize to unity area -----*/
00346   PDF_normalize (p);
00347 
00348 
00349   return;
00350 }
00351 
00352 
00353 /*---------------------------------------------------------------------------*/
00354 /*
00355   Determine the range of values in the input short array.
00356 */
00357 
00358 void PDF_short_range (int npts, short * sarray,
00359                        short * min_val, short * max_val)
00360 {
00361   int ipt;
00362 
00363 
00364   *min_val = sarray[0];
00365   *max_val = sarray[0];
00366 
00367   for (ipt = 1;  ipt < npts;  ipt++)
00368     {
00369       if (sarray[ipt] < *min_val)   *min_val = sarray[ipt];
00370       if (sarray[ipt] > *max_val)   *max_val = sarray[ipt];
00371     }
00372 
00373   return;
00374 }
00375 
00376 
00377 /*---------------------------------------------------------------------------*/
00378 /*
00379   Determine the range of values in the input float array.
00380 */
00381 
00382 void PDF_float_range (int npts, float * farray,
00383                       float * min_val, float * max_val)
00384 {
00385   int ipt;
00386 
00387 
00388   *min_val = farray[0];
00389   *max_val = farray[0];
00390 
00391   for (ipt = 1;  ipt < npts;  ipt++)
00392     {
00393       if (farray[ipt] < *min_val)   *min_val = farray[ipt];
00394       if (farray[ipt] > *max_val)   *max_val = farray[ipt];
00395     }
00396 
00397   return;
00398 }
00399 
00400 
00401 /*---------------------------------------------------------------------------*/
00402 /*
00403   Estimate the pdf corresponding to the input short array.
00404 */
00405 
00406 void PDF_short_to_pdf (int npts, short * sarray, pdf * p)
00407 {
00408   const int MIN_COUNT = 5;
00409   const int MIN_BINS  = 5;
00410   int ipt, ibin, count;
00411   float * fbin = NULL;
00412   int num_bins;
00413   short lower_lim, upper_lim;
00414   char message[80];
00415 
00416 
00417   /*----- Make histogram of input short array -----*/
00418   PDF_short_range (npts, sarray, &lower_lim, &upper_lim);
00419   num_bins = upper_lim - lower_lim + 1;
00420   if (num_bins < MIN_BINS)
00421     {
00422       sprintf (message, "histogram contains only %d bins", num_bins);
00423       PDF_error (message);
00424     }
00425 
00426   fbin = (float *) malloc (sizeof(float) * num_bins);  PDF_MTEST (fbin);
00427   for (ibin = 0;  ibin < num_bins;  ibin++)
00428     fbin[ibin] = 0.0;
00429     
00430   count = 0;
00431   for (ipt = 0;  ipt < npts;  ipt++)
00432     {
00433       ibin = sarray[ipt] - lower_lim;
00434       if ((ibin >= 0) && (ibin < num_bins))
00435         {
00436           fbin[ibin] += 1.0;
00437           count++;
00438         }
00439     }
00440 
00441 
00442   /*----- Check for too few points -----*/
00443   if (count < MIN_COUNT)
00444     {
00445       sprintf (message, "histogram contains only %d points", count);
00446       PDF_error (message);
00447     }
00448 
00449 
00450   /*----- Create PDF -----*/
00451   PDF_create (num_bins, fbin, (float) lower_lim, (float) upper_lim, p);
00452 
00453 
00454   /*----- Release memory -----*/
00455   free (fbin);   fbin = NULL;
00456 
00457 
00458   return;
00459 }
00460  
00461  
00462 /*---------------------------------------------------------------------------*/
00463 /*
00464   Estimate the pdf corresponding to the input float array.
00465 */
00466 void PDF_float_to_pdf (int npts, float * farray, int num_bins, pdf * p)
00467 {
00468   const int MIN_COUNT = 5;
00469   const int MIN_BINS  = 5;
00470   int ipt, ibin, count;
00471   float * fbin = NULL;
00472   float width;
00473   float min_val, max_val;
00474   char message[80];
00475 
00476 
00477   /*----- Make histogram of input float array -----*/
00478   if (num_bins < MIN_BINS)
00479     {
00480       sprintf (message, "histogram contains only %d bins", num_bins);
00481       PDF_error (message);
00482     }
00483 
00484   fbin = (float *) malloc (sizeof(float) * num_bins);   PDF_MTEST (fbin);
00485   for (ibin = 0;  ibin < num_bins;  ibin++)
00486     fbin[ibin] = 0.0;
00487   
00488   PDF_float_range (npts, farray, &min_val, &max_val);
00489   width = (max_val - min_val) / num_bins;
00490  
00491   count = 0;
00492   for (ipt = 0;  ipt < npts;  ipt++)
00493     {
00494       ibin = (farray[ipt] - min_val) / width;
00495       if ((ibin >= 0) && (ibin < num_bins))
00496         {
00497           fbin[ibin] += 1.0;
00498           count++;
00499         }
00500     }
00501 
00502 
00503   /*----- Check for too few points -----*/
00504   if (count < MIN_COUNT)
00505     {
00506       sprintf (message, "histogram contains only %d points", count);
00507       PDF_error (message);
00508     }
00509 
00510 
00511   /*----- Create PDF -----*/
00512   PDF_create (num_bins, fbin, min_val, max_val, p);
00513 
00514 
00515   /*----- Release memory -----*/
00516   free (fbin);   fbin = NULL;
00517 
00518 
00519   return;
00520 }
00521  
00522  
00523 /*---------------------------------------------------------------------------*/
00524 /*
00525   Find extrema of pdf function.
00526 */
00527 
00528 void PDF_find_extrema (pdf p, int * num_min, int * pdf_min, 
00529                        int * num_max, int * pdf_max)
00530 {
00531   int ibin;
00532   int i;
00533 
00534 
00535   *num_min = 0;
00536   *num_max = 0;
00537 
00538   for (ibin = 1;  ibin < p.nbin-1;  ibin++)
00539     {
00540       if ((p.prob[ibin] < p.prob[ibin-1]) && (p.prob[ibin] < p.prob[ibin+1]))
00541         {
00542           pdf_min[*num_min] = ibin;
00543           (*num_min)++;
00544         }
00545 
00546       if ((p.prob[ibin] > p.prob[ibin-1]) && (p.prob[ibin] > p.prob[ibin+1]))
00547         {
00548           pdf_max[*num_max] = ibin;
00549           (*num_max)++;
00550         }
00551     }
00552 
00553   if( !quiet ){
00554    printf ("\nExtrema of PDF: \n");
00555    printf ("\nNum Local Min = %d \n", *num_min);
00556    for (i = 0;  i < *num_min;  i++)
00557      {
00558        ibin = pdf_min[i]; 
00559        printf ("x[%3d] = %8.3f   p[%3d] = %12.6f \n", 
00560                ibin, PDF_ibin_to_xvalue(p, ibin), ibin, p.prob[ibin]);
00561      }
00562 
00563    printf ("\nNum Local Max = %d \n", *num_max);
00564    for (i = 0;  i < *num_max;  i++)
00565      {
00566        ibin = pdf_max[i]; 
00567        printf ("x[%3d] = %8.3f   p[%3d] = %12.6f \n", 
00568                ibin, PDF_ibin_to_xvalue(p, ibin), ibin, p.prob[ibin]);
00569      }
00570   }
00571 
00572 }
00573 
00574 
00575 /*---------------------------------------------------------------------------*/
00576 /*
00577   Find bimodality of pdf function (if possible).
00578 */
00579 
00580 int PDF_find_bimodal (pdf p, int * gmax, int * wmax)
00581 {
00582   const int NPTS = 12;
00583   int * pdf_min = NULL, * pdf_max = NULL;
00584   int num_min, num_max;
00585   int imax, temp;
00586   
00587 
00588   pdf_min = (int *) malloc (sizeof(int) * p.nbin);
00589   pdf_max = (int *) malloc (sizeof(int) * p.nbin);
00590   
00591   PDF_find_extrema (p, &num_min, pdf_min, &num_max, pdf_max);
00592 
00593 
00594   if (num_max >= 2)
00595     {
00596       if (p.prob[pdf_max[1]] >= p.prob[pdf_max[0]])
00597         {
00598           *wmax = pdf_max[1];
00599           *gmax = pdf_max[0];
00600         }
00601       else
00602         {
00603           *wmax = pdf_max[0];
00604           *gmax = pdf_max[1];
00605         }
00606       
00607       if (num_max > 2)
00608         {
00609           for (imax = 2;  imax < num_max;  imax++)
00610             {
00611               if (p.prob[pdf_max[imax]] >= p.prob[*wmax])
00612                 {
00613                   *gmax = *wmax;
00614                   *wmax = pdf_max[imax];
00615                 }
00616               else if (p.prob[pdf_max[imax]] >= p.prob[*gmax])
00617                 {
00618                   *gmax = pdf_max[imax];
00619                 }
00620             }
00621         }
00622 
00623       if (*gmax > *wmax)
00624         {
00625           temp = *gmax;
00626           *gmax = *wmax;
00627           *wmax = temp;
00628         }
00629 
00630     }  /* end if (num_max >= 2) */
00631 
00632 
00633   free (pdf_min);   pdf_min = NULL;
00634   free (pdf_max);   pdf_max = NULL;
00635   
00636 
00637   if (num_max < 2)  return (0);
00638   else              return (1);
00639 
00640 }
00641 
00642 
00643 /*---------------------------------------------------------------------------*/
00644 
00645 
00646 
00647 
00648 
00649