contrib/mul/mbl/mbl_histogram.cxx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_histogram.cxx
00002 #include "mbl_histogram.h"
00003 //:
00004 // \file
00005 // \brief Simple object to build histogram from supplied data.
00006 // \author Tim Cootes
00007 
00008 #include <vcl_cmath.h>
00009 #include <vcl_iostream.h>
00010 #include <vcl_cassert.h>
00011 #include <vsl/vsl_vector_io.h>
00012 
00013 mbl_histogram::mbl_histogram()
00014 {
00015   clear();
00016 }
00017 
00018 // Construct with given number of bins over given range
00019 mbl_histogram::mbl_histogram(double x_lo, double x_hi, int n_bins)
00020 {
00021   set_bins(x_lo,x_hi,n_bins);
00022 }
00023 
00024 //: Define number and size of bins
00025 void mbl_histogram::set_bins(double xlo, double xhi, int n_bins)
00026 {
00027   assert(n_bins>0);
00028   assert(xhi>xlo);
00029 
00030   bins_.resize(n_bins+1);
00031   freq_.resize(n_bins);
00032 
00033   double dx = (xhi-xlo)/n_bins;
00034   for (int i=0;i<=n_bins;++i) bins_[i]=xlo+i*dx;
00035   clear();
00036 }
00037 
00038 void mbl_histogram::clear()
00039 {
00040   n_obs_ = 0;
00041   n_below_ = 0;
00042   n_above_ = 0;
00043   for (unsigned int i=0;i<freq_.size();++i) freq_[i]=0;
00044 }
00045 
00046 void mbl_histogram::obs(double v)
00047 {
00048   n_obs_++;
00049   if (v<bins_[0])
00050   {
00051     n_below_++;
00052     return;
00053   }
00054 
00055   int n = freq_.size();
00056 
00057   for (int i=1;i<=n;++i)
00058   {
00059     if (v<bins_[i])
00060     {
00061       freq_[i-1]++;
00062       return;
00063     }
00064   }
00065 
00066   // Not in any bin
00067   n_above_++;
00068 }
00069 
00070 
00071 const double MAX_ERROR = 1.0e-8;
00072 
00073 //: Test for equality
00074 bool mbl_histogram::operator==(const mbl_histogram& s) const
00075 {
00076   if (s.n_bins()!=n_bins()) return false;
00077   if (s.n_obs_ != n_obs_) return false;
00078   if (s.n_below_!=n_below_) return false;
00079   if (s.n_above_!=n_above_) return false;
00080 
00081   int n = n_bins();
00082   for (int i=0;i<n;++i)
00083     if (s.freq_[i]!=freq_[i]) return false;
00084   for (int i=0;i<=n;++i)
00085     if (vcl_fabs(s.bins_[i]-bins_[i])>MAX_ERROR) return false;
00086 
00087   return true;
00088 }
00089 
00090 //: Version number for I/O
00091 short mbl_histogram::version_no() const
00092 {
00093   return 1;
00094 }
00095 
00096 void mbl_histogram::b_write(vsl_b_ostream& bfs) const
00097 {
00098   vsl_b_write(bfs,version_no());
00099   vsl_b_write(bfs,n_obs_);
00100   vsl_b_write(bfs,n_below_);
00101   vsl_b_write(bfs,n_above_);
00102   vsl_b_write(bfs,bins_);
00103   vsl_b_write(bfs,freq_);
00104 }
00105 
00106 void mbl_histogram::b_read(vsl_b_istream& bfs)
00107 {
00108   if (!bfs) return;
00109 
00110   short file_version_no;
00111   vsl_b_read(bfs,file_version_no);
00112 
00113   switch (file_version_no)
00114   {
00115    case 1:
00116     vsl_b_read(bfs,n_obs_);
00117     vsl_b_read(bfs,n_below_);
00118     vsl_b_read(bfs,n_above_);
00119     vsl_b_read(bfs,bins_);
00120     vsl_b_read(bfs,freq_);
00121     break;
00122    default:
00123     vcl_cerr << "I/O ERROR: mbl_histogram::b_read(vsl_b_istream&)\n"
00124              << "           Unknown version number "<< file_version_no << '\n';
00125     bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00126     return;
00127   }
00128 }
00129 
00130 void mbl_histogram::print_summary(vcl_ostream& os) const
00131 {
00132   os << "mbl_histogram: ";
00133   if (n_bins()==0) { os<< "No bins defined."; return; }
00134 
00135   int n = n_bins();
00136 
00137   if (n_obs_==0)
00138     os << "No samples.";
00139   else
00140   {
00141     os << n_obs_ << " observations.\n"
00142        << "    < "<<bins_[0]<<"   "<<n_below_<<'\n';
00143     for (int i=0;i<n;++i)
00144       os<<"  ["<<bins_[i]<<','<<bins_[i+1]<<")  "<<freq_[i]<<'\n';
00145     os << "   >= "<<bins_[n]<<"   "<<n_above_<<'\n';
00146   }
00147 }
00148 
00149 //: Write out histogram probabilities to a named file
00150 //  Format: (bin-centre) prob     (one per line)
00151 // \return true if successful
00152 bool mbl_histogram::write_probabilities(const char* path)
00153 {
00154   int n = n_bins();
00155   if (n==0) return false;
00156 
00157   vcl_ofstream ofs(path);
00158   if (!ofs) return false;
00159   for (int i=0;i<n_bins();++i)
00160   {
00161     double dx=vcl_fabs(bins_[i+1]-bins_[i]);
00162     ofs<<0.5*(bins_[i]+bins_[i+1])<<"  "<<double(freq_[i])/(dx*n_obs_)<<'\n';
00163   }
00164   ofs.close();
00165   return true;
00166 }
00167 
00168 //: Write out cumulative probability distribution to a named file
00169 //  Format: (bin-centre) sum_prob     (one per line)
00170 // \return true if successful
00171 bool mbl_histogram::write_cdf(const char* path)
00172 {
00173   int n = n_bins();
00174   if (n==0) return false;
00175 
00176   vcl_ofstream ofs(path);
00177   if (!ofs) return false;
00178   int sum=n_below_;
00179   for (int i=0;i<n_bins();++i)
00180   {
00181     sum+=freq_[i];
00182     ofs<<0.5*(bins_[i]+bins_[i+1])<<"  "<<double(sum)/n_obs_<<'\n';
00183   }
00184   ofs.close();
00185   return true;
00186 }
00187 
00188 vcl_ostream& operator<<(vcl_ostream& os, const mbl_histogram& histo)
00189 {
00190   histo.print_summary(os);
00191   return os;
00192 }
00193 
00194   //: Stream output operator for class reference
00195 void vsl_print_summary(vcl_ostream& os,const mbl_histogram& histo)
00196 {
00197   histo.print_summary(os);
00198 }
00199 
00200   //: Binary file stream output operator for class reference
00201 void vsl_b_write(vsl_b_ostream& bfs, const mbl_histogram& h)
00202 {
00203   h.b_write(bfs);
00204 }
00205 
00206   //: Binary file stream input operator for class reference
00207 void vsl_b_read(vsl_b_istream& bfs, mbl_histogram& h)
00208 {
00209   h.b_read(bfs);
00210 }