merger.h

Go to the documentation of this file.

#ifndef MERGER_H
#define MERGER_H


#include "../estim/mixtures.h"

namespace bdm
{
using std::string;

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

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

                MERGER_METHOD METHOD;
                static const MERGER_METHOD DFLT_METHOD;
                
                double beta;
                static const double DFLT_beta;
                
                eEmp eSmp;

                bool DBG;

                it_file* dbg_file;
        public:

                merger_base () : compositepdf() {DBG = false;dbg_file = NULL;}

                merger_base (const Array<mpdf*> &S, bool own=false);

                void set_sources (const Array<mpdf*> &Sources, bool own) {
                        compositepdf::set_elements (Sources,own);
                        Nsources=mpdfs.length();
                        //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) ;
                        };
                }
                void set_support (const Array<vec> &XYZ, const int dimsize) {
                        set_support(XYZ,dimsize*ones_i(XYZ.length()));
                }
                void set_support (const Array<vec> &XYZ, const ivec &gridsize) {
                        int dim = XYZ.length();  //check with internal dim!!
                        Npoints = prod (gridsize); 
                        eSmp.set_parameters (Npoints, false);
                        Array<vec> &samples = eSmp._samples();
                        eSmp._w() = ones (Npoints) / Npoints; //unifrom size of bins
                        //set samples
                        ivec ind = zeros_i (dim);      //indeces of dimensions in for cycle;
                        vec smpi (dim);            // ith sample
                        vec steps =zeros(dim);            // ith sample
                        // first corner
                        for (int j = 0; j < dim; j++) {
                                smpi (j) = XYZ (j) (0); /* beginning of the interval*/ 
                                it_assert(gridsize(j)!=0.0,"Zeros in gridsize!");
                                steps (j) = ( XYZ(j)(1)-smpi(j) )/gridsize(j);
                        }
                        // fill samples
                        for (int i = 0; i < Npoints; i++) {
                                // copy 
                                samples(i) = smpi; 
                                // go through all dimensions
                                for (int j = 0;j < dim;j++) {
                                        if (ind (j) == gridsize (j) - 1) { //j-th index is full
                                                ind (j) = 0; //shift back
                                                smpi(j) = XYZ(j)(0);
                                                
                                                if (i<Npoints-1) {
                                                        ind (j + 1) ++; //increase the next dimension;
                                                        smpi(j+1) += steps(j+1);
                                                        break;
                                                }
                                                
                                        } else {
                                                ind (j) ++; 
                                                smpi(j) +=steps(j);
                                                break;
                                        }
                                }
                        }
                }
                void set_debug_file (const string fname) {
                        if (DBG) delete dbg_file;
                        dbg_file = new it_file (fname);
                        if (dbg_file) DBG = true;
                }
                void set_method (MERGER_METHOD MTH=DFLT_METHOD, double beta0 = DFLT_beta) {
                        METHOD = MTH;
                        beta = beta0;
                }
                void set_support (const epdf &overall, int N) {
                        eSmp.set_statistics (&overall, N);
                        Npoints = N;
                }

                virtual ~merger_base() {
                        for (int i = 0; i < Nsources; i++) {
                                delete dls (i);
                                delete zdls (i);
                        }
                        if (DBG) delete dbg_file;
                };


                virtual void merge () {
                        validate();

                        //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 w_nn=merge_points (lW);
                                vec wtmp = exp (w_nn-max(w_nn));
                                //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 < Npoints; 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 < Npoints; 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 < Nsources; i++) {
                                tmp += w (i) * pow (S (i), 2);
                        }
                        return tmp -pow (mean(), 2);
                }


                eEmp& _Smp() {return eSmp;}

                void from_setting (const Setting& set) {
                        // get support
                        // find which method to use
                        string meth_str;
                        UI::get<string> (meth_str, set, "method", UI::compulsory);
                        if (!strcmp (meth_str.c_str(), "arithmetic"))
                                set_method (ARITHMETIC);
                        else {
                                if (!strcmp (meth_str.c_str(), "geometric"))
                                        set_method (GEOMETRIC);
                                else if (!strcmp (meth_str.c_str(), "lognormal")) {
                                        set_method (LOGNORMAL);
                                        set.lookupValue( "beta",beta);
                                }
                        }
                        string dbg_file;
                        if (UI::get(dbg_file, set, "dbg_file"))
                                set_debug_file(dbg_file);
                        //validate() - not used
                }

                void validate() {
                        it_assert (eSmp._w().length() > 0, "Empty support, use set_support().");
                        it_assert (dim == eSmp._samples() (0).length(), "Support points and rv are not compatible!");
                        it_assert (isnamed(),"mergers must be named");
                }
};
UIREGISTER(merger_base);

class merger_mix : public merger_base
{
        protected:
                MixEF Mix;
                int Ncoms;
                double effss_coef;
                int stop_niter;
                
                static const int DFLT_Ncoms;
                static const double DFLT_effss_coef;

        public:
                merger_mix () {};
                merger_mix (const Array<mpdf*> &S,bool own=false) {set_sources (S,own);};
                void set_sources (const Array<mpdf*> &S, bool own) {
                        merger_base::set_sources (S,own);
                        Nsources = S.length();
                }
                void set_parameters (int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef) {
                        Ncoms = Ncoms0;
                        effss_coef = effss_coef0;
                }


                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);
                }

                MixEF& _Mix() {return Mix;}
                emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
                void from_setting(const Setting& set){
                        merger_base::from_setting(set);
                        set.lookupValue("ncoms",Ncoms);
                        set.lookupValue("effss_coef",effss_coef);
                        set.lookupValue("stop_niter",stop_niter);
                }
                

};
UIREGISTER(merger_mix);

}

#endif // MER_H

---