Go to the documentation of this file.00001
00002 #include "pdf1d_gaussian_builder.h"
00003
00004
00005
00006 #include <vcl_cassert.h>
00007 #include <vcl_string.h>
00008 #include <vcl_cstdlib.h>
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
00017
00018
00019 pdf1d_gaussian_builder::pdf1d_gaussian_builder()
00020 : min_var_(1.0e-6)
00021 {
00022 }
00023
00024
00025
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
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
00053
00054 void pdf1d_gaussian_builder::set_min_var(double min_var)
00055 {
00056 min_var_ = min_var;
00057 }
00058
00059
00060
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
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
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());
00148
00149
00150
00151
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
00173
00174
00175 vcl_string pdf1d_gaussian_builder::is_a() const
00176 {
00177 return vcl_string("pdf1d_gaussian_builder");
00178 }
00179
00180
00181
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
00191
00192
00193 short pdf1d_gaussian_builder::version_no() const
00194 {
00195 return 1;
00196 }
00197
00198
00199
00200
00201
00202 pdf1d_builder* pdf1d_gaussian_builder::clone() const
00203 {
00204 return new pdf1d_gaussian_builder(*this);
00205 }
00206
00207
00208
00209
00210
00211 void pdf1d_gaussian_builder::print_summary(vcl_ostream& os) const
00212 {
00213 os << "Min. var.: "<< min_var_;
00214 }
00215
00216
00217
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
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);
00245 return;
00246 }
00247 }
00248