contrib/mul/pdf1d/pdf1d_gaussian_builder.cxx
Go to the documentation of this file.
00001 // This is mul/pdf1d/pdf1d_gaussian_builder.cxx
00002 #include "pdf1d_gaussian_builder.h"
00003 //:
00004 // \file
00005 
00006 #include <vcl_cassert.h>
00007 #include <vcl_string.h>
00008 #include <vcl_cstdlib.h> // vcl_abort()
00009 
00010 #include <mbl/mbl_data_wrapper.h>
00011 #include <mbl/mbl_data_array_wrapper.h>
00012 #include <pdf1d/pdf1d_gaussian.h>
00013 #include <pdf1d/pdf1d_calc_mean_var.h>
00014 
00015 //=======================================================================
00016 // Dflt ctor
00017 //=======================================================================
00018 
00019 pdf1d_gaussian_builder::pdf1d_gaussian_builder()
00020     : min_var_(1.0e-6)
00021 {
00022 }
00023 
00024 //=======================================================================
00025 // Destructor
00026 //=======================================================================
00027 
00028 pdf1d_gaussian_builder::~pdf1d_gaussian_builder()
00029 {
00030 }
00031 
00032 //=======================================================================
00033 
00034 pdf1d_gaussian& pdf1d_gaussian_builder::gaussian(pdf1d_pdf& model) const
00035 {
00036   // require a pdf1d_gaussian
00037   assert(model.is_class("pdf1d_gaussian"));
00038   return static_cast<pdf1d_gaussian&>(model);
00039 }
00040 
00041 pdf1d_pdf* pdf1d_gaussian_builder::new_model() const
00042 {
00043   return new pdf1d_gaussian;
00044 }
00045 
00046 vcl_string pdf1d_gaussian_builder::new_model_type() const
00047 {
00048   return vcl_string("pdf1d_gaussian");
00049 }
00050 
00051 //=======================================================================
00052 //: Define lower threshold on variance for built models
00053 //=======================================================================
00054 void pdf1d_gaussian_builder::set_min_var(double min_var)
00055 {
00056   min_var_ = min_var;
00057 }
00058 
00059 //=======================================================================
00060 //: Get lower threshold on variance for built models
00061 //=======================================================================
00062 double pdf1d_gaussian_builder::min_var() const
00063 {
00064   return min_var_;
00065 }
00066 
00067 void pdf1d_gaussian_builder::build(pdf1d_pdf& model, double mean) const
00068 {
00069   pdf1d_gaussian& g = gaussian(model);
00070   g.set(mean,min_var_);
00071 }
00072 
00073 //: Build gaussian from n elements in data[i]
00074 void pdf1d_gaussian_builder::build_from_array(pdf1d_pdf& model, const double* data, int n) const
00075 {
00076   if (n<2)
00077   {
00078     vcl_cerr<<"pdf1d_gaussian_builder::build_from_array() Too few samples available.\n";
00079     vcl_abort();
00080   }
00081 
00082   double m,v;
00083   pdf1d_calc_mean_var(m,v,data,n);
00084 
00085   if (v<min_var_) v=min_var_;
00086 
00087   pdf1d_gaussian& g = gaussian(model);
00088   g.set(m,v);
00089 }
00090 
00091 void pdf1d_gaussian_builder::build(pdf1d_pdf& model, mbl_data_wrapper<double>& data) const
00092 {
00093   pdf1d_gaussian& g = gaussian(model);
00094 
00095   int n_samples = data.size();
00096 
00097   if (n_samples<2)
00098   {
00099     vcl_cerr<<"pdf1d_gaussian_builder::build() Too few samples available.\n";
00100     vcl_abort();
00101   }
00102 
00103   if (data.is_class("mbl_data_array_wrapper<T>"))
00104   {
00105     // Use more efficient build_from_array algorithm
00106     mbl_data_array_wrapper<double>& data_array =
00107                      static_cast<mbl_data_array_wrapper<double>&>(data);
00108     build_from_array(model,data_array.data(),n_samples);
00109     return;
00110   }
00111 
00112   double sum = 0;
00113   double sum_sq = 0;
00114 
00115   data.reset();
00116   for (int i=0; i++<n_samples; data.next())
00117   {
00118     double x = data.current();
00119     sum += x;
00120     sum_sq += x*x;
00121   }
00122 
00123   double m = sum/n_samples;
00124   double v = (sum_sq - m*sum)/(n_samples-1);
00125   if (v<min_var_) v=min_var_;
00126 
00127   g.set(m,v);
00128 }
00129 
00130 void pdf1d_gaussian_builder::weighted_build(pdf1d_pdf& model,
00131                                             mbl_data_wrapper<double>& data,
00132                                             const vcl_vector<double>& wts) const
00133 {
00134   pdf1d_gaussian& g = gaussian(model);
00135 
00136   int n_samples = data.size();
00137 
00138   if (n_samples<2)
00139   {
00140     vcl_cerr<<"pdf1d_gaussian_builder::build() Too few samples available.\n";
00141     vcl_abort();
00142   }
00143 
00144   double sum = 0;
00145   double w_sum = 0;
00146   double sum_sq = 0;
00147   const double* w = &(wts.front()); // cannot use wts.begin() since that's an iterator, not a pointer
00148 
00149 
00150   // Inefficient to go through twice.
00151   // Fix this one day.
00152 
00153   data.reset();
00154   for (int i=0;i<n_samples;++i,data.next())
00155   {
00156     double x = data.current();
00157     double wi = w[i];
00158     sum += wi*x;
00159     sum_sq += wi*x*x;
00160     w_sum += wi;
00161   }
00162 
00163   double m = sum/w_sum;
00164   double v = sum_sq/w_sum -m*m;
00165   v *= n_samples/double(n_samples-1);
00166   if (v<min_var_) v=min_var_;
00167 
00168   g.set(m,v);
00169 }
00170 
00171 //=======================================================================
00172 // Method: is_a
00173 //=======================================================================
00174 
00175 vcl_string pdf1d_gaussian_builder::is_a() const
00176 {
00177   return vcl_string("pdf1d_gaussian_builder");
00178 }
00179 
00180 //=======================================================================
00181 // Method: is_class
00182 //=======================================================================
00183 
00184 bool pdf1d_gaussian_builder::is_class(vcl_string const& s) const
00185 {
00186   return pdf1d_builder::is_class(s) || s==pdf1d_gaussian_builder::is_a();
00187 }
00188 
00189 //=======================================================================
00190 // Method: version_no
00191 //=======================================================================
00192 
00193 short pdf1d_gaussian_builder::version_no() const
00194 {
00195   return 1;
00196 }
00197 
00198 //=======================================================================
00199 // Method: clone
00200 //=======================================================================
00201 
00202 pdf1d_builder* pdf1d_gaussian_builder::clone() const
00203 {
00204   return new pdf1d_gaussian_builder(*this);
00205 }
00206 
00207 //=======================================================================
00208 // Method: print
00209 //=======================================================================
00210 
00211 void pdf1d_gaussian_builder::print_summary(vcl_ostream& os) const
00212 {
00213   os << "Min. var.: "<< min_var_;
00214 }
00215 
00216 //=======================================================================
00217 // Method: save
00218 //=======================================================================
00219 
00220 void pdf1d_gaussian_builder::b_write(vsl_b_ostream& bfs) const
00221 {
00222   vsl_b_write(bfs,version_no());
00223   vsl_b_write(bfs,min_var_);
00224 }
00225 
00226 //=======================================================================
00227 // Method: load
00228 //=======================================================================
00229 
00230 void pdf1d_gaussian_builder::b_read(vsl_b_istream& bfs)
00231 {
00232   if (!bfs) return;
00233 
00234   short version;
00235   vsl_b_read(bfs,version);
00236   switch (version)
00237   {
00238     case (1):
00239       vsl_b_read(bfs,min_var_);
00240       break;
00241     default:
00242       vcl_cerr << "I/O ERROR: vsl_b_read(vsl_b_istream&, pdf1d_gaussian_builder &)\n"
00243                << "           Unknown version number "<< version << '\n';
00244       bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00245       return;
00246   }
00247 }
00248