root/library/bdm/stat/merger.h @ 1085

/*!
  \file
  \brief Mergers for combination of pdfs
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef MERGER_H
#define MERGER_H


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

namespace bdm {
using std::string;

//!Merging methods
enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3};

/*!  \brief Abstract common class for mergers

It defines only common interface of access to sources and operation merge(); The result of merging is available via function merger()
*/
class MergerBase: public root{
        public:         
                //! weights of the sources -- no restrictions
                vec w;
                //! source pdfs -- conditioning allowed
                Array<shared_ptr<epdf> > sources;
                //! merge 
                virtual void merge() NOT_IMPLEMENTED_VOID;
                //! 
                //! check if all epdfs are on the same support
                virtual void set_sources(const Array<shared_ptr<pdf> > &A) {
                        int n=A.length();
                        
                        bdm_assert(n>0,"merger has no sources to merge");
                        // check if weights are ok
                        if (w.length()!=n) {
                                w = ones(n)/n;
                        }
                        
                        // check compatibility of sources -- no types are needed
                        int dim0 = A(0)->dimension();
                        for (int i=0; i<n; i++){
                                const epdf *si=dynamic_cast<const epdf*>(A(i).get());
                                if (si){
                                        sources(i) = dynamic_cast<epdf*>(A(i)->_copy());
                                }
                                bdm_assert(sources(i)->dimension()==dim0, "Merger: Incompatible dimensions of sources");
                                if (sources(i)->isnamed()){
                                        //check if rv match
                                }
                        }
                }
                //! return the result of merging
                virtual const epdf& merger()=0;
};

//! \brief Merger into normal density, max entropy approximator for 2 moments (mean+variance)
template<class sq_T>
class ENormMerger: public MergerBase{
        protected:
                //!internal epdf
                enorm<sq_T> iepdf;
        public:
                ENormMerger():method(GEOMETRIC){};
                MERGER_METHOD method;
                void merge(){
                        int n = sources.length();
                        int dim = sources(0)->dimension();
                        sq_T Var(zeros(dim));
                        vec mea=zeros(dim);
                        // go through all sources
                        epdf * si;
                        for (int i=0; i<n; i++){
                                si = dynamic_cast<epdf*>(sources(i).get());  // transform pdf into epdf

                                if (!si) {bdm_error("Can't merger this type of pdf: " + sources(i)->to_string());}
                                
                                sq_T Ci = si->covariance();
                                sq_T wCi = Ci; wCi*=w(i);
                                vec mi = si->mean();
                                switch (method) {
                                        case ARITHMETIC:{
                                                Var += wCi;
                                                Var.opupdt(mi,w(i)); // + mean * mean'
                                                mea += w(i)*mi;
                                                break;
                                        }
                                        case GEOMETRIC:{
                                                Var += wCi;
                                                sq_T iCi;
                                                Ci.inv(iCi);
                                                mea += iCi.to_mat()*w(i)*mi;
                                                break;
                                        }
                                        default:
                                                bdm_error("Method not implemneted");
                                }
                        }
                        // post pocessing
                        switch (method) {
                                case ARITHMETIC:
                                        iepdf._R()=Var;
                                        iepdf._R().opupdt(mea,-1.0); // Var - mean * mean'
                                        iepdf._mu()=mea;
                                        break;
                                case GEOMETRIC:
                                        iepdf._R()=Var;
                                        iepdf._mu() = Var.to_mat()*mea; // mean=sqrt(Var)*mea
                                        break;
                                default:
                                        bdm_error("Method not implemneted");
                        }
                }
                //! access function
                const epdf& merger(){return iepdf;}
};

/*!
@brief Base class (interface) for for general combination of pdfs on discrete support

Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.

The merged pdfs are expected to be of the form:
 \f[ f(x_i|y_i),  i=1..n \f]
where the resulting merger is a density on \f$ \cup [x_i,y_i] \f$ .
Note that all variables will be joined.

As a result of this feature, each source must be extended to common support
\f[ f(z_i|y_i,x_i) f(x_i|y_i) f(y_i)  i=1..n \f]
where \f$ z_i \f$ accumulate variables that were not in the original source.
These extensions are calculated on-the-fly.

However, these operations can not be performed in general. Hence, this class merges only sources on common support, \f$ y_i={}, z_i={}, \forall i \f$.
For merging of more general cases, use offsprings merger_mix and merger_grid.

\TODO use sources for extensions
*/
class MergerDiscrete : public MergerBase {
protected:
        //! partial sources
        Array<shared_ptr< pdf > > part_sources;
        
    //! Data link for each pdf in pdfs
    Array<datalink_m2e*> dls;

    //! Array of rvs that are not modelled by pdfs at all, \f$ z_i \f$
    Array<RV> rvzs;

    //! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
    Array<datalink_m2e*> zdls;

    //! number of support points
    int Npoints;

    //! number of sources
    int Nsources;

    //! switch of the methoh used for merging
    MERGER_METHOD METHOD;
    //! Default for METHOD
    static const MERGER_METHOD DFLT_METHOD;

    //!Prior on the log-normal merging model
    double beta;
    //! default for beta
    static const double DFLT_beta;

    //! Projection to empirical density (could also be piece-wise linear)
    eEmp eSmp;

    //! debug or not debug
    bool DBG;

    //! debugging file
    it_file* dbg_file;
public:
    //! \name Constructors
    //! @{

    //! Default constructor
    MergerDiscrete () : Npoints ( 0 ), Nsources ( 0 ), DBG ( false ), dbg_file ( 0 ) {
    }

    //!Constructor from sources
    MergerDiscrete ( const Array<shared_ptr<pdf> > &S );

    //! Function setting the main internal structures
    void set_sources ( const Array<shared_ptr<pdf> > &PartSources );

    //! Set support points from rectangular grid
    void set_support ( rectangular_support &Sup );

    //! Set support points from dicrete grid
    void set_support ( discrete_support &Sup ) {
        Npoints = Sup.points();
        eSmp.set_parameters ( Sup._Spoints() );
        eSmp.validate();
    }


    //! set debug file
    void set_debug_file ( const string fname ) {
        if ( DBG ) delete dbg_file;
        dbg_file = new it_file ( fname );
        DBG = ( dbg_file != 0 );
    }

    //! Set internal parameters used in approximation
    void set_method ( MERGER_METHOD MTH = GEOMETRIC, double beta0 = 1.2 ) {
        METHOD = MTH;
        beta = beta0;
    }
    //! Set support points from a pdf by drawing N samples
    void set_support ( const epdf &overall, int N ) {
        eSmp.set_statistics ( overall, N );
        Npoints = N;
    }

    //! Destructor
    virtual ~MergerDiscrete() {
        for ( int i = 0; i < Nsources; i++ ) {
            delete dls ( i );
            delete zdls ( i );
        }
        if ( DBG ) delete dbg_file;
    };
    //!@}

    //! \name Mathematical operations
    //!@{

    //!Merge given sources in given points
    void merge ();

    //! Merge log-likelihood values in points using method specified by parameter METHOD
    vec merge_points ( mat &lW );

  //!@}

    //! \name Access to attributes
    //! @{

    //! Access function
    eEmp& _Smp() {
        return eSmp;
    }

        eEmp & merger() {return eSmp;}

    /*! Create object from the following structure

    \code
    class = 'merger_base';
    method = 'arithmetic' or 'geometric' or 'lognormal';   % name of the model used for merging 
    --- compulsory only for lognormal merging model ---
    beta = x;                                              % scalar prior on the log-normal merging model
    --- optional fields ---
    dbg_file = '...';                                      % name of the file to store debug informations
    --- inherited fields ---
    bdm::epdf::from_setting
    \endcode
    */
    void from_setting ( const Setting& set );

    void to_setting  (Setting  &set) const ;

    void validate() ;
    //!@}
};
UIREGISTER ( MergerDiscrete );
SHAREDPTR ( MergerDiscrete );

//! \brief Merger using importance sampling with mixture proposal density
class merger_mix : public MergerDiscrete {
protected:
    //!Internal mixture of EF models
    MixEF Mix;
    //!Number of components in a mixture
    int Ncoms;
    //! coefficient of resampling [0,1]
    double effss_coef;
    //! stop after niter iterations
    int stop_niter;

    //! default value for Ncoms
    static const int DFLT_Ncoms;
    //! default value for efss_coef;
    static const double DFLT_effss_coef;

public:
    //!\name Constructors
    //!@{
    merger_mix () : Ncoms ( 0 ), effss_coef ( 0 ), stop_niter ( 0 ) { }

    merger_mix ( const Array<shared_ptr<pdf> > &S ) :
        Ncoms ( 0 ), effss_coef ( 0 ), stop_niter ( 0 ) {
        set_sources ( S );
    }

    //! Set sources and prepare all internal structures
    void set_sources ( const Array<shared_ptr<pdf> > &S ) {
        MergerDiscrete::set_sources ( S );
        //Nsources = S.length();
    }

    //! Set internal parameters used in approximation
    void set_parameters ( int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef ) {
        Ncoms = Ncoms0;
        effss_coef = effss_coef0;
    }
    //!@}

    //! \name Mathematical operations
    //!@{

    //!Merge values using mixture approximation
    void merge ();

    //! sample from the approximating mixture
    vec sample () const {
        return Mix.posterior().sample();
    }
    //! loglikelihood computed on mixture models
    double evallog ( const vec &yt ) const {
        vec dtf = ones ( yt.length() + 1 );
        dtf.set_subvector ( 0, yt );
        return Mix.logpred ( dtf, vec(0) );
    }
    //!@}

    //!\name Access functions
    //!@{
//! Access function
    MixEF& _Mix() {
        return Mix;
    }
    //! Access function
    emix* proposal() {
        emix* tmp = Mix.epredictor();
        tmp->set_rv ( merger()._rv() );
        return tmp;
    }
    /*! Create object from the following structure

    \code
    class = 'merger_mix';
    --- optional fields ---
    ncoms = [];        % number of components in a mixture
    effss_coef = [];   % coefficient of resampling from interval [0,1]
    stop_niter = [];   % number of interations
    --- inherited fields ---
    bdm::merger_base::from_setting
    \endcode    
    If the optional fields are not given, they will be filled as follows:
    \code
    ncoms = 10;       
    effss_coef = 0.9;  
    stop_niter = 10;  
    \endcode

    */
    void from_setting ( const Setting& set );

    void to_setting  (Setting  &set) const;
    void validate();

//      const emix &merger(){return *Mpred;}
    //! @}
};
UIREGISTER ( merger_mix );
SHAREDPTR ( merger_mix );

}

#endif // MER_H