Go to the documentation of this file.00001
00002 #include "mbl_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_histogram::mbl_histogram()
00014 {
00015 clear();
00016 }
00017
00018
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
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
00067 n_above_++;
00068 }
00069
00070
00071 const double MAX_ERROR = 1.0e-8;
00072
00073
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
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);
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
00150
00151
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
00169
00170
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
00195 void vsl_print_summary(vcl_ostream& os,const mbl_histogram& histo)
00196 {
00197 histo.print_summary(os);
00198 }
00199
00200
00201 void vsl_b_write(vsl_b_ostream& bfs, const mbl_histogram& h)
00202 {
00203 h.b_write(bfs);
00204 }
00205
00206
00207 void vsl_b_read(vsl_b_istream& bfs, mbl_histogram& h)
00208 {
00209 h.b_read(bfs);
00210 }