contrib/mul/mbl/mbl_wt_histogram.cxx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_wt_histogram.cxx
00002 #include "mbl_wt_histogram.h"
00003 //:
00004 // \file
00005 // \brief Simple object to build histogram from supplied data, with weights
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_wt_histogram::mbl_wt_histogram()
00014 {
00015   clear();
00016 }
00017 
00018 // Construct with given number of bins over given range
00019 mbl_wt_histogram::mbl_wt_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_wt_histogram::set_bins(double xlo, double xhi, int n_bins)
00026 {
00027   assert(n_bins>0);
00028   assert(xhi>xlo);
00029 
00030   wt_sum_.resize(n_bins);
00031 
00032   dx_ = (xhi-xlo)/n_bins;
00033   xlo_ = xlo;
00034   clear();
00035 }
00036 
00037 void mbl_wt_histogram::clear()
00038 {
00039   n_obs_ = 0;
00040   total_wt_=0.0;
00041   wt_below_ = 0;
00042   wt_above_ = 0;
00043   for (unsigned int i=0;i<wt_sum_.size();++i) wt_sum_[i]=0;
00044 }
00045 
00046 void mbl_wt_histogram::obs(double v, double wt)
00047 {
00048   n_obs_++;
00049   total_wt_ += wt;
00050   if (v<xlo_)
00051   {
00052     wt_below_+=wt;
00053     return;
00054   }
00055 
00056   // v-xlo_ >= 0
00057   unsigned int i = (unsigned int)((v-xlo_)/dx_);
00058 
00059   if (i<wt_sum_.size()) wt_sum_[i]+=wt;
00060   else                  wt_above_+=wt;
00061 }
00062 
00063 
00064 const double MAX_ERROR = 1.0e-8;
00065 
00066 //: Test for equality
00067 bool mbl_wt_histogram::operator==(const mbl_wt_histogram& s) const
00068 {
00069   if (s.n_bins()!=n_bins()) return false;
00070   if (s.n_obs_ != n_obs_) return false;
00071   if (vcl_fabs(s.wt_below_-wt_below_)>MAX_ERROR) return false;
00072   if (vcl_fabs(s.wt_above_-wt_above_)>MAX_ERROR) return false;
00073   if (vcl_fabs(s.xlo_-xlo_)>MAX_ERROR) return false;
00074   if (vcl_fabs(s.dx_-dx_)>MAX_ERROR) return false;
00075 
00076   int n = n_bins();
00077   for (int i=0;i<n;++i)
00078     if (vcl_fabs(s.wt_sum_[i]-wt_sum_[i])>MAX_ERROR) return false;
00079 
00080   return true;
00081 }
00082 
00083 //: Version number for I/O
00084 short mbl_wt_histogram::version_no() const
00085 {
00086   return 1;
00087 }
00088 
00089 void mbl_wt_histogram::b_write(vsl_b_ostream& bfs) const
00090 {
00091   vsl_b_write(bfs,version_no());
00092   vsl_b_write(bfs,n_obs_);
00093   vsl_b_write(bfs,total_wt_);
00094   vsl_b_write(bfs,wt_below_);
00095   vsl_b_write(bfs,wt_above_);
00096   vsl_b_write(bfs,xlo_);
00097   vsl_b_write(bfs,dx_);
00098   vsl_b_write(bfs,wt_sum_);
00099 }
00100 
00101 void mbl_wt_histogram::b_read(vsl_b_istream& bfs)
00102 {
00103   if (!bfs) return;
00104 
00105   short file_version_no;
00106   vsl_b_read(bfs,file_version_no);
00107 
00108   switch (file_version_no)
00109   {
00110    case 1:
00111     vsl_b_read(bfs,n_obs_);
00112     vsl_b_read(bfs,total_wt_);
00113     vsl_b_read(bfs,wt_below_);
00114     vsl_b_read(bfs,wt_above_);
00115     vsl_b_read(bfs,xlo_);
00116     vsl_b_read(bfs,dx_);
00117     vsl_b_read(bfs,wt_sum_);
00118     break;
00119    default:
00120     vcl_cerr << "I/O ERROR: mbl_wt_histogram::b_read(vsl_b_istream&)\n"
00121              << "           Unknown version number "<< file_version_no << '\n';
00122     bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00123     return;
00124   }
00125 }
00126 
00127 void mbl_wt_histogram::print_summary(vcl_ostream& os) const
00128 {
00129   os << "mbl_wt_histogram: ";
00130   if (n_bins()==0) { os<< "No bins defined."; return; }
00131 
00132   int n = n_bins();
00133 
00134   if (n_obs_==0)
00135     os << "No samples.";
00136   else
00137   {
00138     os << n_obs_ << " observations.\n"
00139        << "    < "<<xlo_<<"   "<<wt_below_<<'\n';
00140     for (int i=0;i<n;++i)
00141       os<<"  ["<<xlo_+i*dx_<<','<<xlo_+(i+1)*dx_<<")  "<<wt_sum_[i]<<'\n';
00142     os << "   >= "<<xlo_+n*dx_<<"   "<<wt_above_<<'\n';
00143   }
00144 }
00145 
00146 //: Write out histogram probabilities to a named file
00147 //  Format: (bin-centre) prob     (one per line)
00148 // \return true if successful
00149 bool mbl_wt_histogram::write_probabilities(const char* path)
00150 {
00151   int n = n_bins();
00152   if (n==0) return false;
00153 
00154   vcl_ofstream ofs(path);
00155   if (!ofs) return false;
00156   for (int i=0;i<n;++i)
00157   {
00158     ofs<<xlo_+(i+0.5)*dx_<<"  "<<double(wt_sum_[i])/total_wt_<<'\n';
00159   }
00160   ofs.close();
00161   return true;
00162 }
00163 
00164 vcl_ostream& operator<<(vcl_ostream& os, const mbl_wt_histogram& histo)
00165 {
00166   histo.print_summary(os);
00167   return os;
00168 }
00169 
00170   //: Stream output operator for class reference
00171 void vsl_print_summary(vcl_ostream& os,const mbl_wt_histogram& histo)
00172 {
00173   histo.print_summary(os);
00174 }
00175 
00176   //: Binary file stream output operator for class reference
00177 void vsl_b_write(vsl_b_ostream& bfs, const mbl_wt_histogram& h)
00178 {
00179   h.b_write(bfs);
00180 }
00181 
00182   //: Binary file stream input operator for class reference
00183 void vsl_b_read(vsl_b_istream& bfs, mbl_wt_histogram& h)
00184 {
00185   h.b_read(bfs);
00186 }