- Main Page
- Related Pages
- Modules
- Classes
- Files

- File List
- File Members

# merger.h

Go to the documentation of this file.

    00001
    00013 #ifndef MERGER_H
    00014 #define MERGER_H
    00015
    00016
    00017 #include "../estim/mixtures.h"
    00018
    00019 namespace bdm
    00020 {
    00021 using std::string;
    00022
    00024 enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3};
    00025
    00045 class merger_base : public compositepdf, public epdf
    00046 {
    00047         protected:
    00049                 Array<datalink_m2e*> dls;
    00051                 Array<RV> rvzs;
    00053                 Array<datalink_m2e*> zdls;
    00055                 int Npoints;
    00057                 int Nsources;
    00058
    00060                 MERGER_METHOD METHOD;
    00062                 static const MERGER_METHOD DFLT_METHOD;
    00063
    00065                 double beta;
    00067                 static const double DFLT_beta;
    00068
    00070                 eEmp eSmp;
    00071
    00073                 bool DBG;
    00074
    00076                 it_file* dbg_file;
    00077         public:
    00080
    00082                 merger_base () : compositepdf() {DBG = false;dbg_file = NULL;};
    00084                 merger_base (const Array<mpdf*> &S, bool own=false) {set_sources (S,own);};
    00086                 void set_sources (const Array<mpdf*> &Sources, bool own) {
    00087                         compositepdf::set_elements (Sources,own);
    00088                         Nsources=mpdfs.length();
    00089                         //set sizes
    00090                         dls.set_size (Sources.length());
    00091                         rvzs.set_size (Sources.length());
    00092                         zdls.set_size (Sources.length());
    00093
    00094                         rv = getrv (/* checkoverlap = */ false);
    00095                         RV rvc; setrvc (rv, rvc);  // Extend rv by rvc!
    00096                         // join rv and rvc - see descriprion
    00097                         rv.add (rvc);
    00098                         // get dimension
    00099                         dim = rv._dsize();
    00100
    00101                         // create links between sources and common rv
    00102                         RV xytmp;
    00103                         for (int i = 0;i < mpdfs.length();i++) {
    00104                                 //Establich connection between mpdfs and merger
    00105                                 dls (i) = new datalink_m2e;
    00106                                 dls (i)->set_connection (mpdfs (i)->_rv(), mpdfs (i)->_rvc(), rv);
    00107
    00108                                 // find out what is missing in each mpdf
    00109                                 xytmp = mpdfs (i)->_rv();
    00110                                 xytmp.add (mpdfs (i)->_rvc());
    00111                                 // z_i = common_rv-xy
    00112                                 rvzs (i) = rv.subt (xytmp);
    00113                                 //establish connection between extension (z_i|x,y)s and common rv
    00114                                 zdls (i) = new datalink_m2e; zdls (i)->set_connection (rvzs (i), xytmp, rv) ;
    00115                         };
    00116                 }
    00118                 void set_support (const Array<vec> &XYZ, const int dimsize) {
    00119                         set_support(XYZ,dimsize*ones_i(XYZ.length()));
    00120                 }
    00122                 void set_support (const Array<vec> &XYZ, const ivec &gridsize) {
    00123                         int dim = XYZ.length();  //check with internal dim!!
    00124                         Npoints = prod (gridsize);
    00125                         eSmp.set_parameters (Npoints, false);
    00126                         Array<vec> &samples = eSmp._samples();
    00127                         eSmp._w() = ones (Npoints) / Npoints; //unifrom size of bins
    00128                         //set samples
    00129                         ivec ind = zeros_i (dim);      //indeces of dimensions in for cycle;
    00130                         vec smpi (dim);            // ith sample
    00131                         vec steps =zeros(dim);            // ith sample
    00132                         // first corner
    00133                         for (int j = 0; j < dim; j++) {
    00134                                 smpi (j) = XYZ (j) (0); /* beginning of the interval*/
    00135                                 it_assert(gridsize(j)!=0.0,"Zeros in gridsize!");
    00136                                 steps (j) = ( XYZ(j)(1)-smpi(j) )/gridsize(j);
    00137                         }
    00138                         // fill samples
    00139                         for (int i = 0; i < Npoints; i++) {
    00140                                 // copy
    00141                                 samples(i) = smpi;
    00142                                 // go through all dimensions
    00143                                 for (int j = 0;j < dim;j++) {
    00144                                         if (ind (j) == gridsize (j) - 1) { //j-th index is full
    00145                                                 ind (j) = 0; //shift back
    00146                                                 smpi(j) = XYZ(j)(0);
    00147
    00148                                                 if (i<Npoints-1) {
    00149                                                         ind (j + 1) ++; //increase the next dimension;
    00150                                                         smpi(j+1) += steps(j+1);
    00151                                                         break;
    00152                                                 }
    00153
    00154                                         } else {
    00155                                                 ind (j) ++;
    00156                                                 smpi(j) +=steps(j);
    00157                                                 break;
    00158                                         }
    00159                                 }
    00160                         }
    00161                 }
    00163                 void set_debug_file (const string fname) {
    00164                         if (DBG) delete dbg_file;
    00165                         dbg_file = new it_file (fname);
    00166                         if (dbg_file) DBG = true;
    00167                 }
    00169                 void set_method (MERGER_METHOD MTH=DFLT_METHOD, double beta0 = DFLT_beta) {
    00170                         METHOD = MTH;
    00171                         beta = beta0;
    00172                 }
    00174                 void set_support (const epdf &overall, int N) {
    00175                         eSmp.set_statistics (&overall, N);
    00176                         Npoints = N;
    00177                 }
    00178
    00180                 virtual ~merger_base() {
    00181                         for (int i = 0; i < Nsources; i++) {
    00182                                 delete dls (i);
    00183                                 delete zdls (i);
    00184                         }
    00185                         if (DBG) delete dbg_file;
    00186                 };
    00188
    00191
    00193                 virtual void merge () {
    00194                         validate();
    00195
    00196                         //check if sources overlap:
    00197                         bool OK = true;
    00198                         for (int i = 0;i < mpdfs.length(); i++) {
    00199                                 OK &= (rvzs (i)._dsize() == 0); // z_i is empty
    00200                                 OK &= (mpdfs (i)->_rvc()._dsize() == 0); // y_i is empty
    00201                         }
    00202
    00203                         if (OK) {
    00204                                 mat lW = zeros (mpdfs.length(), eSmp._w().length());
    00205
    00206                                 vec emptyvec (0);
    00207                                 for (int i = 0; i < mpdfs.length(); i++) {
    00208                                         for (int j = 0; j < eSmp._w().length(); j++) {
    00209                                                 lW (i, j) = mpdfs (i)->evallogcond (eSmp._samples() (j), emptyvec);
    00210                                         }
    00211                                 }
    00212
    00213                                 vec w_nn=merge_points (lW);
    00214                                 vec wtmp = exp (w_nn-max(w_nn));
    00215                                 //renormalize
    00216                                 eSmp._w() = wtmp / sum (wtmp);
    00217                         } else {
    00218                                 it_error ("Sources are not compatible - use merger_mix");
    00219                         }
    00220                 };
    00221
    00222
    00224                 vec merge_points (mat &lW);
    00225
    00226
    00229                 vec mean() const {
    00230                         const Vec<double> &w = eSmp._w();
    00231                         const Array<vec> &S = eSmp._samples();
    00232                         vec tmp = zeros (dim);
    00233                         for (int i = 0; i < Npoints; i++) {
    00234                                 tmp += w (i) * S (i);
    00235                         }
    00236                         return tmp;
    00237                 }
    00238                 mat covariance() const {
    00239                         const vec &w = eSmp._w();
    00240                         const Array<vec> &S = eSmp._samples();
    00241
    00242                         vec mea = mean();
    00243
    00244                         cout << sum (w) << "," << w*w << endl;
    00245
    00246                         mat Tmp = zeros (dim, dim);
    00247                         for (int i = 0; i < Npoints; i++) {
    00248                                 Tmp += w (i) * outer_product (S (i), S (i));
    00249                         }
    00250                         return Tmp -outer_product (mea, mea);
    00251                 }
    00252                 vec variance() const {
    00253                         const vec &w = eSmp._w();
    00254                         const Array<vec> &S = eSmp._samples();
    00255
    00256                         vec tmp = zeros (dim);
    00257                         for (int i = 0; i < Nsources; i++) {
    00258                                 tmp += w (i) * pow (S (i), 2);
    00259                         }
    00260                         return tmp -pow (mean(), 2);
    00261                 }
    00263
    00266
    00268                 eEmp& _Smp() {return eSmp;}
    00269
    00271                 void from_setting (const Setting& set) {
    00272                         // get support
    00273                         // find which method to use
    00274                         string meth_str;
    00275                         UI::get<string> (meth_str, set, "method");
    00276                         if (!strcmp (meth_str.c_str(), "arithmetic"))
    00277                                 set_method (ARITHMETIC);
    00278                         else {
    00279                                 if (!strcmp (meth_str.c_str(), "geometric"))
    00280                                         set_method (GEOMETRIC);
    00281                                 else if (!strcmp (meth_str.c_str(), "lognormal")) {
    00282                                         set_method (LOGNORMAL);
    00283                                         set.lookupValue( "beta",beta);
    00284                                 }
    00285                         }
    00286                         if (set.exists("dbg_file")){
    00287                                 string dbg_file;
    00288                                 UI::get<string> (dbg_file, set, "dbg_file");
    00289                                 set_debug_file(dbg_file);
    00290                         }
    00291                         //validate() - not used
    00292                 }
    00293
    00294                 void validate() {
    00295                         it_assert (eSmp._w().length() > 0, "Empty support, use set_support().");
    00296                         it_assert (dim == eSmp._samples() (0).length(), "Support points and rv are not compatible!");
    00297                         it_assert (isnamed(),"mergers must be named");
    00298                 }
    00300 };
    00301 UIREGISTER(merger_base);
    00302
    00303 class merger_mix : public merger_base
    00304 {
    00305         protected:
    00307                 MixEF Mix;
    00309                 int Ncoms;
    00311                 double effss_coef;
    00313                 int stop_niter;
    00314
    00316                 static const int DFLT_Ncoms;
    00318                 static const double DFLT_effss_coef;
    00319
    00320         public:
    00323                 merger_mix () {};
    00324                 merger_mix (const Array<mpdf*> &S,bool own=false) {set_sources (S,own);};
    00326                 void set_sources (const Array<mpdf*> &S, bool own) {
    00327                         merger_base::set_sources (S,own);
    00328                         Nsources = S.length();
    00329                 }
    00331                 void set_parameters (int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef) {
    00332                         Ncoms = Ncoms0;
    00333                         effss_coef = effss_coef0;
    00334                 }
    00336
    00339
    00341                 void merge ();
    00342
    00344                 vec sample () const { return Mix.posterior().sample();}
    00346                 double evallog (const vec &dt) const {
    00347                         vec dtf = ones (dt.length() + 1);
    00348                         dtf.set_subvector (0, dt);
    00349                         return Mix.logpred (dtf);
    00350                 }
    00352
    00356                 MixEF& _Mix() {return Mix;}
    00358                 emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
    00360                 void from_setting(const Setting& set){
    00361                         merger_base::from_setting(set);
    00362                         set.lookupValue("ncoms",Ncoms);
    00363                         set.lookupValue("effss_coef",effss_coef);
    00364                         set.lookupValue("stop_niter",stop_niter);
    00365                 }
    00366
    00368
    00369 };
    00370 UIREGISTER(merger_mix);
    00371
    00372 }
    00373
    00374 #endif // MER_H