merger.h

Go to the documentation of this file.

#ifndef MER_H
#define MER_H


#include "mixef.h"

namespace bdm
{
using std::string;

enum MERGER_METHOD {ARITHMETIC = 0, GEOMETRIC = 1, LOGNORMAL = 2};

class merger_base : public compositepdf, public epdf
{
        protected:
                Array<datalink_m2e*> dls;
                Array<RV> rvzs;
                Array<datalink_m2e*> zdls;
                int Ns;

                MERGER_METHOD METHOD;
                double beta;

                eEmp eSmp;

                bool DBG;

                it_file* dbg;
        public:

                merger_base () : compositepdf() {};
                merger_base (const Array<mpdf*> &S) {set_sources (S);};
                void set_sources (const Array<mpdf*> &Sources) {
                        compositepdf::set_elements (Sources);
                        //set sizes
                        dls.set_size (Sources.length());
                        rvzs.set_size (Sources.length());
                        zdls.set_size (Sources.length());

                        rv = getrv (/* checkoverlap = */ false);
                        RV rvc; setrvc (rv, rvc);  // Extend rv by rvc!
                        // join rv and rvc - see descriprion
                        rv.add (rvc);
                        // get dimension
                        dim = rv._dsize();

                        // create links between sources and common rv
                        RV xytmp;
                        for (int i = 0;i < mpdfs.length();i++) {
                                //Establich connection between mpdfs and merger
                                dls (i) = new datalink_m2e;
                                dls (i)->set_connection (mpdfs (i)->_rv(), mpdfs (i)->_rvc(), rv);

                                // find out what is missing in each mpdf
                                xytmp = mpdfs (i)->_rv();
                                xytmp.add (mpdfs (i)->_rvc());
                                // z_i = common_rv-xy
                                rvzs (i) = rv.subt (xytmp);
                                //establish connection between extension (z_i|x,y)s and common rv
                                zdls (i) = new datalink_m2e; zdls (i)->set_connection (rvzs (i), xytmp, rv) ;
                        };
                        //Set Default values of parameters
                        DBG = false;
                        METHOD = LOGNORMAL;
                        beta = 2.0;
                }
                void set_debug_file (const string fname) { if (DBG) delete dbg; dbg = new it_file (fname); if (dbg) DBG = true;}
                void set_parameters (double beta0, MERGER_METHOD MTH) {beta = beta0;METHOD = MTH;}
                void set_support(epdf &overall, int N){eSmp.set_statistics(&overall,N);Ns=N;}
                
                virtual ~merger_base() {
                        for (int i = 0; i < mpdfs.length(); i++) {
                                delete dls (i);
                                delete zdls (i);
                        }
                        if (DBG) delete dbg;
                };
                
                
                void merge () {
                        if (eSmp._w().length() ==0) {it_error("Empty support points use set_support" );}
                        //check if sources overlap:
                        bool OK = true;
                        for (int i = 0;i < mpdfs.length(); i++) {
                                OK &= (rvzs (i)._dsize() == 0); // z_i is empty
                                OK &= (mpdfs (i)->_rvc()._dsize() == 0); // y_i is empty
                        }

                        if (OK) {
                                mat lW = zeros (mpdfs.length(), eSmp._w().length());

                                vec emptyvec (0);
                                for (int i = 0; i < mpdfs.length(); i++) {
                                        for (int j = 0; j < eSmp._w().length(); j++) {
                                                lW (i, j) = mpdfs (i)->evallogcond (eSmp._samples() (j), emptyvec);
                                        }
                                }

                                vec wtmp = exp (merge_points (lW));
                                //renormalize
                                eSmp._w() = wtmp / sum (wtmp);
                        } else {
                                it_error("Sources are not compatible - use merger_mix");
                        }
                        ;
                };


                vec merge_points (mat &lW);
                
                
                vec mean() const {
                        const Vec<double> &w = eSmp._w();
                        const Array<vec> &S = eSmp._samples();
                        vec tmp = zeros (dim);
                        for (int i = 0; i < Ns; i++) {
                                tmp += w (i) * S (i);
                        }
                        return tmp;
                }
                mat covariance() const {
                        const vec &w = eSmp._w();
                        const Array<vec> &S = eSmp._samples();

                        vec mea = mean();

                        cout << sum (w) << "," << w*w << endl;

                        mat Tmp = zeros (dim, dim);
                        for (int i = 0; i < Ns; i++) {
                                Tmp += w (i) * outer_product (S (i), S (i));
                        }
                        return Tmp -outer_product (mea, mea);
                }
                vec variance() const {
                        const vec &w = eSmp._w();
                        const Array<vec> &S = eSmp._samples();

                        vec tmp = zeros (dim);
                        for (int i = 0; i < Ns; i++) {
                                tmp += w (i) * pow (S (i), 2);
                        }
                        return tmp -pow (mean(), 2);
                }


                eEmp& _Smp() {return eSmp;}
                
};

class merger_mix : public merger_base
{
        protected:
                MixEF Mix;
                int Nc;
                double effss_coef;

        public:
                merger_mix () {};
                merger_mix (const Array<mpdf*> &S) {set_sources(S);};
                void set_sources (const Array<mpdf*> &S) {
                        merger_base::set_sources(S);
                        //Set Default values of parameters
                        Ns = 100;
                        Nc = 10;
                        Mix.set_method (EM);
                }
                void set_parameters (double beta0, int Ns0, int Nc0, double effss_coef0 = 0.5) {
                        beta = beta0;
                        Ns = Ns0;
                        Nc = Nc0;
                        effss_coef = effss_coef0;
                        eSmp.set_parameters (Ns0, false);
                }
                void init() { 
                        Array<vec> Smps (Ns);
                        //Gibbs sampling
                        for (int i = 0;i < mpdfs.length();i++) {Smps (i) = zeros (0);}
                }
                void merge ();

                vec sample () const { return Mix.posterior().sample();}
                double evallog (const vec &dt) const {
                        vec dtf = ones (dt.length() + 1);
                        dtf.set_subvector (0, dt);
                        return Mix.logpred (dtf);
                }
                vec mean() const {
                        const Vec<double> &w = eSmp._w();
                        const Array<vec> &S = eSmp._samples();
                        vec tmp = zeros (dim);
                        for (int i = 0; i < Ns; i++) {
                                tmp += w (i) * S (i);
                        }
                        return tmp;
                }
                mat covariance() const {
                        const vec &w = eSmp._w();
                        const Array<vec> &S = eSmp._samples();

                        vec mea = mean();

                        cout << sum (w) << "," << w*w << endl;

                        mat Tmp = zeros (dim, dim);
                        for (int i = 0; i < Ns; i++) {
                                Tmp += w (i) * outer_product (S (i), S (i));
                        }
                        return Tmp -outer_product (mea, mea);
                }
                vec variance() const {
                        const vec &w = eSmp._w();
                        const Array<vec> &S = eSmp._samples();

                        vec tmp = zeros (dim);
                        for (int i = 0; i < Ns; i++) {
                                tmp += w (i) * pow (S (i), 2);
                        }
                        return tmp -pow (mean(), 2);
                }

                MixEF& _Mix() {return Mix;}
                emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
                eEmp& _Smp() {return eSmp;}
};

}

#endif // MER_H

---