Go to the documentation of this file.00001
00002 #include "mbl_wt_histogram.h"
00003
00004
00005
00006
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
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
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
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
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
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);
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
00147
00148
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
00171 void vsl_print_summary(vcl_ostream& os,const mbl_wt_histogram& histo)
00172 {
00173 histo.print_summary(os);
00174 }
00175
00176
00177 void vsl_b_write(vsl_b_ostream& bfs, const mbl_wt_histogram& h)
00178 {
00179 h.b_write(bfs);
00180 }
00181
00182
00183 void vsl_b_read(vsl_b_istream& bfs, mbl_wt_histogram& h)
00184 {
00185 h.b_read(bfs);
00186 }