/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
/*---------------------------------------------------------------------------*/
/*
This file contains general purpose routines for input, manipulation, and
output of probability density functions.
File: pdf.c
Author: B. Douglas Ward
Date: 28 January 2000
*/
#ifndef USE_QUIET
# define quiet 0
#endif
/*---------------------------------------------------------------------------*/
/*
Routine to print an error message and stop.
*/
void PDF_error (char * message)
{
printf ("PDF error: %s \n", message);
exit (1);
}
/*---------------------------------------------------------------------------*/
/*
Initialize pdf data structure.
*/
void PDF_initialize (pdf * p)
{
p->nbin = 0;
p->prob = NULL;
p->lower_bnd = 0.0;
p->upper_bnd = 0.0;
p->width = 0.0;
return;
}
/*---------------------------------------------------------------------------*/
/*
Destroy pdf data structure by deallocating memory.
*/
void PDF_destroy (pdf * p)
{
if (p->prob != NULL) free (p->prob);
PDF_initialize (p);
return;
}
/*---------------------------------------------------------------------------*/
/*
Normalize pdf to have unity area.
*/
void PDF_normalize (pdf * p)
{
int ibin;
double sum;
sum = 0.0;
for (ibin = 0; ibin < p->nbin; ibin++)
sum += p->prob[ibin];
for (ibin = 0; ibin < p->nbin; ibin++)
p->prob[ibin] /= sum;
return;
}
/*---------------------------------------------------------------------------*/
/*
Create pdf data structure by allocating memory and initializing values.
*/
void PDF_create (int nbin, float * prob, float lower_bnd, float upper_bnd,
pdf * p)
{
int ibin;
PDF_destroy (p);
p->nbin = nbin;
p->prob = (float *) malloc (sizeof(float) * nbin);
PDF_MTEST (p->prob);
for (ibin = 0; ibin < nbin; ibin++)
p->prob[ibin] = prob[ibin];
p->lower_bnd = lower_bnd;
p->upper_bnd = upper_bnd;
p->width = (upper_bnd - lower_bnd) / (nbin-1);
PDF_normalize (p);
return;
}
/*---------------------------------------------------------------------------*/
/*
Copy pdf data structure from p to pc.
*/
void PDF_copy (pdf p, pdf * pc)
{
PDF_create (p.nbin, p.prob, p.lower_bnd, p.upper_bnd, pc);
return;
}
/*---------------------------------------------------------------------------*/
/*
Convert bin number to value of independent variable at center of the bin.
*/
float PDF_ibin_to_xvalue (pdf p, int ibin)
{
float xvalue;
xvalue = p.lower_bnd + ibin*p.width;
return (xvalue);
}
/*---------------------------------------------------------------------------*/
/*
Convert value of independent variable to bin number.
*/
int PDF_xvalue_to_ibin (pdf p, float xvalue)
{
int ibin;
ibin = floor ( (xvalue - p.lower_bnd) / p.width + 0.5);
return (ibin);
}
/*---------------------------------------------------------------------------*/
/*
Convert value of independent variable to probability.
*/
float PDF_xvalue_to_pvalue (pdf p, float xvalue)
{
int ibin;
float pvalue;
ibin = PDF_xvalue_to_ibin (p, xvalue);
if ((ibin < 0) || (ibin >= p.nbin))
pvalue = 0.0;
else
pvalue = p.prob[ibin];
return (pvalue);
}
/*---------------------------------------------------------------------------*/
/*
Print contents of pdf p to screen.
*/
void PDF_print (pdf p)
{
int ibin;
if( !quiet ){
printf ("Number of bins = %d \n", p.nbin);
printf ("Lower bound = %f \n", p.lower_bnd);
printf ("Upper bound = %f \n", p.upper_bnd);
printf ("Bin width = %f \n", p.width);
/*
printf ("%3s %10.6s %10.6s \n", "i", "x[i]", "p[i]");
for (ibin = 0; ibin < p.nbin; ibin++)
printf ("%3d %10.6f %10.6f \n",
ibin, PDF_ibin_to_xvalue(p, ibin), p.prob[ibin]);
*/
}
return;
}
/*---------------------------------------------------------------------------*/
/*
Print contents of string str and pdf p to screen.
*/
void PDF_sprint (char * str, pdf p)
{
if( quiet ) return ;
printf ("%s \n", str);
PDF_print (p);
return;
}
/*---------------------------------------------------------------------------*/
/*
Write contents of pdf p to specified file.
*/
void PDF_write_file (char * filename, pdf p)
{
int ibin;
FILE * outfile = NULL;
outfile = fopen (filename, "w");
if (!outfile) {
fprintf (stderr, "\n"
"*****************************\n"
"Error:\n"
"Failed to open %s for output.\n"
"Check for write permissions.\n"
"*****************************\n"
"\n", filename);
return;
}
for (ibin = 0; ibin < p.nbin; ibin++)
fprintf (outfile, "%d %f %f \n",
ibin, PDF_ibin_to_xvalue(p, ibin), p.prob[ibin]);
fclose (outfile);
return;
}
/*---------------------------------------------------------------------------*/
/*
Simple smoothing of the pdf estimate.
*/
void PDF_smooth (pdf * p)
{
float * sprob;
int ibin;
sprob = (float *) malloc (sizeof(float) * p->nbin);
sprob[0] = 0.5*(p->prob[0] + p->prob[1]);
sprob[p->nbin-1] = 0.5*(p->prob[p->nbin-2] + p->prob[p->nbin-1]);
for (ibin = 1; ibin < p->nbin-1; ibin++)
sprob[ibin] = 0.25 * (p->prob[ibin-1] + 2*p->prob[ibin] + p->prob[ibin+1]);
free (p->prob);
p->prob = sprob;
PDF_normalize (p);
return;
}
/*---------------------------------------------------------------------------*/
/*
Trim the pdf by removing the extreme lower and extreme upper values.
*/
void PDF_trim (float lower_per, float upper_per, pdf * p)
{
int ibin;
float * fbin = NULL;
float cum_prob;
float lower_bnd, upper_bnd;
int lo_bin, hi_bin;
/*----- Trim lower values -----*/
cum_prob = 0.0;
for (ibin = 0; ibin < p->nbin; ibin++)
{
cum_prob += p->prob[ibin];
p->prob[ibin] = 0.0;
if (cum_prob > lower_per)
{
lo_bin = ibin + 1;
break;
}
}
/*----- Trim upper values -----*/
cum_prob = 0.0;
for (ibin = p->nbin-1; ibin >= 0; ibin--)
{
cum_prob += p->prob[ibin];
p->prob[ibin] = 0.0;
if (cum_prob > 1.0 - upper_per)
{
hi_bin = ibin - 1;
break;
}
}
/*----- Reset lower and upper bounds -----*/
lower_bnd = PDF_ibin_to_xvalue (*p, lo_bin);
upper_bnd = PDF_ibin_to_xvalue (*p, hi_bin);
p->lower_bnd = lower_bnd;
p->upper_bnd = upper_bnd;
p->nbin = hi_bin - lo_bin + 1;
/*----- Copy data -----*/
fbin = (float *) malloc (sizeof(float) * p->nbin);
for (ibin = 0; ibin < p->nbin; ibin++)
fbin[ibin] = p->prob[ibin+lo_bin];
/*----- Reset pointer to data -----*/
free (p->prob);
p->prob = fbin;
/*----- Normalize to unity area -----*/
PDF_normalize (p);
return;
}
/*---------------------------------------------------------------------------*/
/*
Determine the range of values in the input short array.
*/
void PDF_short_range (int npts, short * sarray,
short * min_val, short * max_val)
{
int ipt;
*min_val = sarray[0];
*max_val = sarray[0];
for (ipt = 1; ipt < npts; ipt++)
{
if (sarray[ipt] < *min_val) *min_val = sarray[ipt];
if (sarray[ipt] > *max_val) *max_val = sarray[ipt];
}
return;
}
/*---------------------------------------------------------------------------*/
/*
Determine the range of values in the input float array.
*/
void PDF_float_range (int npts, float * farray,
float * min_val, float * max_val)
{
int ipt;
*min_val = farray[0];
*max_val = farray[0];
for (ipt = 1; ipt < npts; ipt++)
{
if (farray[ipt] < *min_val) *min_val = farray[ipt];
if (farray[ipt] > *max_val) *max_val = farray[ipt];
}
return;
}
/*---------------------------------------------------------------------------*/
/*
Estimate the pdf corresponding to the input short array.
*/
void PDF_short_to_pdf (int npts, short * sarray, pdf * p)
{
const int MIN_COUNT = 5;
const int MIN_BINS = 5;
int ipt, ibin, count;
float * fbin = NULL;
int num_bins;
short lower_lim, upper_lim;
char message[80];
/*----- Make histogram of input short array -----*/
PDF_short_range (npts, sarray, &lower_lim, &upper_lim);
num_bins = upper_lim - lower_lim + 1;
if (num_bins < MIN_BINS)
{
sprintf (message, "histogram contains only %d bins", num_bins);
PDF_error (message);
}
fbin = (float *) malloc (sizeof(float) * num_bins); PDF_MTEST (fbin);
for (ibin = 0; ibin < num_bins; ibin++)
fbin[ibin] = 0.0;
count = 0;
for (ipt = 0; ipt < npts; ipt++)
{
ibin = sarray[ipt] - lower_lim;
if ((ibin >= 0) && (ibin < num_bins))
{
fbin[ibin] += 1.0;
count++;
}
}
/*----- Check for too few points -----*/
if (count < MIN_COUNT)
{
sprintf (message, "histogram contains only %d points", count);
PDF_error (message);
}
/*----- Create PDF -----*/
PDF_create (num_bins, fbin, (float) lower_lim, (float) upper_lim, p);
/*----- Release memory -----*/
free (fbin); fbin = NULL;
return;
}
/*---------------------------------------------------------------------------*/
/*
Estimate the pdf corresponding to the input float array.
*/
void PDF_float_to_pdf (int npts, float * farray, int num_bins, pdf * p)
{
const int MIN_COUNT = 5;
const int MIN_BINS = 5;
int ipt, ibin, count;
float * fbin = NULL;
float width;
float min_val, max_val;
char message[80];
/*----- Make histogram of input float array -----*/
if (num_bins < MIN_BINS)
{
sprintf (message, "histogram contains only %d bins", num_bins);
PDF_error (message);
}
fbin = (float *) malloc (sizeof(float) * num_bins); PDF_MTEST (fbin);
for (ibin = 0; ibin < num_bins; ibin++)
fbin[ibin] = 0.0;
PDF_float_range (npts, farray, &min_val, &max_val);
width = (max_val - min_val) / num_bins;
count = 0;
for (ipt = 0; ipt < npts; ipt++)
{
ibin = (farray[ipt] - min_val) / width;
if ((ibin >= 0) && (ibin < num_bins))
{
fbin[ibin] += 1.0;
count++;
}
}
/*----- Check for too few points -----*/
if (count < MIN_COUNT)
{
sprintf (message, "histogram contains only %d points", count);
PDF_error (message);
}
/*----- Create PDF -----*/
PDF_create (num_bins, fbin, min_val, max_val, p);
/*----- Release memory -----*/
free (fbin); fbin = NULL;
return;
}
/*---------------------------------------------------------------------------*/
/*
Find extrema of pdf function.
*/
void PDF_find_extrema (pdf p, int * num_min, int * pdf_min,
int * num_max, int * pdf_max)
{
int ibin;
int i;
*num_min = 0;
*num_max = 0;
for (ibin = 1; ibin < p.nbin-1; ibin++)
{
if ((p.prob[ibin] < p.prob[ibin-1]) && (p.prob[ibin] < p.prob[ibin+1]))
{
pdf_min[*num_min] = ibin;
(*num_min)++;
}
if ((p.prob[ibin] > p.prob[ibin-1]) && (p.prob[ibin] > p.prob[ibin+1]))
{
pdf_max[*num_max] = ibin;
(*num_max)++;
}
}
if( !quiet ){
printf ("\nExtrema of PDF: \n");
printf ("\nNum Local Min = %d \n", *num_min);
for (i = 0; i < *num_min; i++)
{
ibin = pdf_min[i];
printf ("x[%3d] = %8.3f p[%3d] = %12.6f \n",
ibin, PDF_ibin_to_xvalue(p, ibin), ibin, p.prob[ibin]);
}
printf ("\nNum Local Max = %d \n", *num_max);
for (i = 0; i < *num_max; i++)
{
ibin = pdf_max[i];
printf ("x[%3d] = %8.3f p[%3d] = %12.6f \n",
ibin, PDF_ibin_to_xvalue(p, ibin), ibin, p.prob[ibin]);
}
}
}
/*---------------------------------------------------------------------------*/
/*
Find bimodality of pdf function (if possible).
*/
int PDF_find_bimodal (pdf p, int * gmax, int * wmax)
{
const int NPTS = 12;
int * pdf_min = NULL, * pdf_max = NULL;
int num_min, num_max;
int imax, temp;
pdf_min = (int *) malloc (sizeof(int) * p.nbin);
pdf_max = (int *) malloc (sizeof(int) * p.nbin);
PDF_find_extrema (p, &num_min, pdf_min, &num_max, pdf_max);
if (num_max >= 2)
{
if (p.prob[pdf_max[1]] >= p.prob[pdf_max[0]])
{
*wmax = pdf_max[1];
*gmax = pdf_max[0];
}
else
{
*wmax = pdf_max[0];
*gmax = pdf_max[1];
}
if (num_max > 2)
{
for (imax = 2; imax < num_max; imax++)
{
if (p.prob[pdf_max[imax]] >= p.prob[*wmax])
{
*gmax = *wmax;
*wmax = pdf_max[imax];
}
else if (p.prob[pdf_max[imax]] >= p.prob[*gmax])
{
*gmax = pdf_max[imax];
}
}
}
if (*gmax > *wmax)
{
temp = *gmax;
*gmax = *wmax;
*wmax = temp;
}
} /* end if (num_max >= 2) */
free (pdf_min); pdf_min = NULL;
free (pdf_max); pdf_max = NULL;
if (num_max < 2) return (0);
else return (1);
}
/*---------------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1