root/library/bdm/stat/discrete.h @ 563

/*!
  \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
{

        //! Rectangular support
        //! Support ponits are located inbetween ranges! For example:
        //! For ranges=[0,1] and gridsizes=[1] the support point is 0.5
        class rectangular_support: public root {
                protected:
                        //! Array of boundaries (2D vectors: [begining,end]) for each dimension
                        Array<vec> ranges;
                        //! Number of support points in each dimension
                        ivec gridsizes;
                        //! dimension
                        int dim;
                        //! Number of data points
                        int Npoints;
                        //! active vector for first_vec and next_vec
                        vec actvec;
                        //! indeces of active vector
                        vec actvec_ind;
                        //! length of steps in each dimension
                        vec steps;
                public:
                        //! default constructor
                        rectangular_support() : dim(0), Npoints(0) {
                        }
                        
                        //! set parameters
                        void set_parameters(const Array<vec> &ranges0, const ivec &gridsize0){
                                ranges=ranges0;
                                gridsizes=gridsize0;
                                initialize();
                        }
                        //! Internal functio to set temporaries correctly
                        void initialize() {
                                dim = ranges.length();
                                it_assert_debug(gridsizes.length()==dim,"Incompatible dimensions of input");
                                Npoints = prod(gridsizes);
                                it_assert_debug(Npoints>0,"Wrong input parameters");
                                
                                //precompute steps
                                steps.set_size(dim);
                                for ( int j = 0; j < dim; j++ ) {
                                        steps ( j ) = ( ranges ( j ) ( 1 ) - ranges(j)(0) ) / gridsizes ( j );
                                }
                                actvec.set_size(dim);
                                actvec_ind.set_size(dim);
                        }
                        //! return vector at position given by vector of indeces
                        vec get_vec(const ivec &inds){
                                vec v;
                                for ( int j = 0; j < dim; j++ ) {
                                        it_assert_debug(inds(j)<gridsizes(j), "Index out of bounds");
                                        v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j);
                                }
                                return v;
                        }

                        //! convert dimension indeces into linear index, the indexing is in the same way as in \c next_vec()
                        long linear_index (const ivec inds){
                                long ind=0;
                                it_assert_debug(inds.length()==dim,"Improper indeces inds");
                                int dim_skips=1; // skips in active dimension, in the first dimension, the skips are 1 index per value
                                                
                                for (int j=0; j<dim; j++){ 
                                        ind += dim_skips*(inds(j)); // add shift in linear index caused by this dimension
                                        dim_skips*=gridsizes(j);  // indeces in the next dimension are repeated with period gridsizes(j) times greater that in this dimesion
                                }
                                return ind;
                        } 
                        //! set the first corner to actvec
                        const vec& first_vec(){
                                for ( int j = 0; j < dim; j++ ) {
                                        actvec ( j ) = ranges(j)(0) + 0.5*steps(j);
                                        actvec_ind(j) = 0;
                                }
                                return actvec;
                        }
                        //! Get next active vector, call ONLY after first_vector()!
                        const vec& next_vec() {
                                // go through all dimensions
                                int j=0;
                                while (j<dim) {
                                        if ( actvec_ind ( j ) == gridsizes ( j ) - 1 ) { //j-th index is full
                                                actvec_ind ( j ) = 0; //shift back
                                                actvec ( j ) = ranges ( j ) ( 0 ) + 0.5*steps(j);
                                                j++;
                                        } else {
                                                actvec_ind ( j ) ++;
                                                actvec ( j ) += steps ( j );
                                                break;
                                        }
                                }
                                return actvec;
                        } 
                        
                        ivec nearest_point(const vec &val){
                                ivec inds;
                                inds.set_size(dim);
                                for(int j=0; j<dim; j++){
                                        if (val(j)<ranges(j)(0))
                                                inds(j) = 0;
                                        else {
                                                if (val(j)>ranges(j)(1))
                                                        inds(j) = gridsizes(j)-1;
                                                else {
                                                        inds(j) = ::round(val(j) - ranges(j)(0)/ steps(j));
                                                }
                                        }
                                }
                                return inds;
                        }

                        //! Access function
                        int points() const {return Npoints;}
                        
                        //! access function
                        const vec& _steps() const {return steps;}
                        
                        void from_setting (const Setting &set) {
                                UI::get (ranges , set, "ranges", UI::compulsory);
                                UI::get (gridsizes, set, "gridsizes", UI::compulsory);
                                initialize();
                        }
        };
        UIREGISTER(rectangular_support);
        
        class grid_fnc: public fnc{
                protected:
                        rectangular_support sup;
                        vec values;
                public:
                        //! constructor function
                        void set_support(rectangular_support &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&)){
                                if (sup.points()>0){
                                        values(0) = (*evalptr)(sup.first_vec());
                                        for (int j=1; j<sup.points(); j++){ 
                                                values(j)=(*evalptr)(sup.next_vec());
                                        }
                                }
                        } 
                        //! get value nearest to the given point
                        double nearest_val(const vec &val){return values( sup.linear_index(sup.nearest_point(val)));}
                        
                        vec eval(const vec &val){return vec_1(nearest_val(val));}
        };
        UIREGISTER(grid_fnc);

        //! 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:
                        rectangular_support sup;
                        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