root/library/bdm/stat/discrete.h

/*!
  \file
  \brief Probability distributions for discrete support densities
  \author Vaclav Smidl.

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

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

#ifndef DISCR_H
#define DISCR_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm {

        class support_base: public root {
        protected:
                //! dimension
                int dim;
                //! Number of data points
                int Npoints;
                //! active vector for first_vec and next_vec
                vec actvec;
                //! volume of the active point
                double actvolume;
        public:
                //! default constructor
                support_base() : dim ( 0 ), Npoints ( 0 ) {
                }
                                                
                //! set the first vector to corner and store result in actvec
                virtual const vec& first_vec() NOT_IMPLEMENTED(empty_vec);
                
                //! Get next active vector, call ONLY after first_vector()!
                virtual const vec& next_vec() NOT_IMPLEMENTED(empty_vec);;
                
                //! return active vector, call ONLY after first_vector()!
                virtual const vec& act_vec() {
                        return actvec;
                };
                //! return active vector, call ONLY after first_vector()!
                virtual const double& act_volume() {
                        return actvolume;
                };
                
                //! Access function
                int points() const {
                        return Npoints;
                }
                                
        };
        UIREGISTER ( support_base );
        
//! Rectangular support
//! Support points are located inbetween ranges! For example:
//! For ranges=[0,1] and gridsizes=[1] the support point is 0.5
class rectangular_support: public support_base {
protected:
    //! Array of boundaries (2D vectors: [begining,end]) for each dimension
    Array<vec> ranges;
    //! Number of support points in each dimension
    ivec gridsizes;
        //! indices of active vector
        vec actvec_ind;
        //
        vec steps;
public:
    //! default constructor
    rectangular_support() : support_base() {
    }

    //! set parameters
    void set_parameters ( const Array<vec> &ranges0, const ivec &gridsize0 );

    //! Internal functio to set temporaries correctly
    void initialize();

    //! return vector at position given by vector of indices
    vec get_vec ( const ivec &inds );

    //! convert dimension indices into linear index, the indexing is in the same way as in \c next_vec()
    long linear_index ( const ivec inds );

    //! set the first vector to corner and store result in actvec
    const vec& first_vec();

    //! Get next active vector, call ONLY after first_vector()!
    const vec& next_vec();

    //! return active vector, call ONLY after first_vector()!
    const vec& act_vec() {
        return actvec;
    };

    //! \todo to je asi navic .. v predkovi!
    ivec nearest_point ( const vec &val );

    //! Access function
    int points() const {
        return Npoints;
    }

    //! access function
    //! \todo opet pouze do potomka..
    const vec& _steps() const {
        return steps;
    }

    void from_setting ( const Setting &set );
};
UIREGISTER ( rectangular_support );

//! Discrete support with stored support points
class discrete_support: public support_base {
protected:
    //! storage of support points
    Array<vec> Spoints;
    //! index in iterators
    int idx;
        //! volumes for each point
        vec volumes;
public:
    //! Default constructor
    discrete_support() : Spoints ( 0 ), idx ( 0 ), volumes(0) {}
    //! Access function
    int points() const {
        return Spoints.length();
    }
    //! set the first vector to corner and store result in actvec
    const vec& first_vec() {

                bdm_assert_debug ( Spoints.length() > 0, "Empty support" );
        idx = 0;
                actvolume=volumes(0);
        return Spoints ( idx );
    }
    //! set next vector after calling first_vec()
    const vec& next_vec() {
        bdm_assert_debug ( Spoints.length() > idx - 1, "Out of support points" );
                idx++;
                actvolume=volumes(idx);
                actvec = Spoints ( idx );
        return actvec;
    }

    /*!
    \code
      class = "discrete_support";
      points = ( [1,2..], [2,2..], ...); // list of points
       === OR ===
      epdf = {class="epdf_offspring",...}; // epdf rfom which to sample
      npoints = 100;                     // number of samples
    \endcode
    */
    void from_setting ( const Setting &set );

        void validate(){
                bdm_assert(volumes.length()==Spoints.length(),"dimension mismatch");
                Npoints=Spoints.length();
                if (Npoints>0)
                        dim = Spoints(0).length();
        }
    //! access function
    Array<vec> & _Spoints() {
        return Spoints;
    }
};
UIREGISTER ( discrete_support );

//! Function defined by values on a fixed grid and interpolated inbetween them
class grid_fnc: public fnc {
protected:
    //! grid - function support
    shared_ptr<support_base> sup;
    //! function values on the grid
    vec values;
public:
    //! constructor function
    void set_support ( shared_ptr<support_base> &sup0 ) {
        sup = sup0;
        values = zeros ( sup->points() );
    }
    //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
    void set_values ( double ( *evalptr ) ( const vec& ) );

    //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
    void set_values ( const epdf &ep );

    //! get value nearest to the given point
    double nearest_val ( const vec &val ) {
                NOT_IMPLEMENTED(0.0);
//        return values ( sup.linear_index ( sup.nearest_point ( val ) ) );
    }

        double integrate() ;

    vec eval ( const vec &val ) {
        return vec_1 ( nearest_val ( val ) );
    }
    const vec & _values() const {
        return values;
    }
};
UIREGISTER ( grid_fnc );

//! \brief Piecewise constant pdf on rectangular support
//!
//! Each point on the grid represents a centroid around which the density is constant.
//! This is a trivial point-mass density where all points have the same mass.
class egrid: public epdf {
protected:
    //! support of the pdf - grid
    rectangular_support sup;
    //! values at the grid
    vec values;
public:
    //! we assume that evallog is not called too often otherwise we should cache log(values)
    double evallog ( const vec &val ) {
        return log ( values ( sup.linear_index ( sup.nearest_point ( val ) ) ) );
    }
};
}
#endif //DISCR_H