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

/*!
  \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(){};
                        
                        //! 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);
                        }

                        //! 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
                                for ( int j = 0; j < dim; j++ ) {
                                        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);

                                                //! 
                                                if ( j+1<dim)
                                                        if (actvec_ind ( j + 1 )<gridsizes(j+1)-1){
                                                                actvec_ind ( j + 1 ) ++; //increase the next dimension;
                                                                actvec ( j + 1 ) += steps ( j + 1 );
                                                                break;
                                                        }

                                        } else {
                                                actvec_ind ( j ) ++;
                                                actvec ( j ) += steps ( j );
                                                break;
                                        }
                                }
                                return actvec;
                        } 
                        //! 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);
}
#endif //DISCR_H