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

Revision 569, 6.8 kB (checked in by smidl, 15 years ago)

new object discrete_support, merger adapted to accept this input as well

Line 
1/*!
2  \file
3  \brief Probability distributions for discrete support densities
4  \author Vaclav Smidl.
5
6  -----------------------------------
7  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9  Using IT++ for numerical operations
10  -----------------------------------
11*/
12
13#ifndef DISCR_H
14#define DISCR_H
15
16
17#include "../shared_ptr.h"
18#include "../base/bdmbase.h"
19#include "../math/chmat.h"
20
21namespace bdm
22{
23
24        //! Rectangular support
25        //! Support points are located inbetween ranges! For example:
26        //! For ranges=[0,1] and gridsizes=[1] the support point is 0.5
27        class rectangular_support: public root {
28                protected:
29                        //! Array of boundaries (2D vectors: [begining,end]) for each dimension
30                        Array<vec> ranges;
31                        //! Number of support points in each dimension
32                        ivec gridsizes;
33                        //! dimension
34                        int dim;
35                        //! Number of data points
36                        int Npoints;
37                        //! active vector for first_vec and next_vec
38                        vec actvec;
39                        //! indeces of active vector
40                        vec actvec_ind;
41                        //! length of steps in each dimension
42                        vec steps;
43                public:
44                        //! default constructor
45                        rectangular_support() : dim(0), Npoints(0) {
46                        }
47                       
48                        //! set parameters
49                        void set_parameters(const Array<vec> &ranges0, const ivec &gridsize0){
50                                ranges=ranges0;
51                                gridsizes=gridsize0;
52                                initialize();
53                        }
54                        //! Internal functio to set temporaries correctly
55                        void initialize() {
56                                dim = ranges.length();
57                                bdm_assert_debug(gridsizes.length() == dim, "Incompatible dimensions of input");
58                                Npoints = prod(gridsizes);
59                                bdm_assert_debug(Npoints > 0, "Wrong input parameters");
60                               
61                                //precompute steps
62                                steps.set_size(dim);
63                                for ( int j = 0; j < dim; j++ ) {
64                                        steps ( j ) = ( ranges ( j ) ( 1 ) - ranges(j)(0) ) / gridsizes ( j );
65                                }
66                                actvec.set_size(dim);
67                                actvec_ind.set_size(dim);
68                        }
69                        //! return vector at position given by vector of indeces
70                        vec get_vec(const ivec &inds){
71                                vec v ( dim );
72                                for ( int j = 0; j < dim; j++ ) {
73                                        bdm_assert_debug(inds(j) < gridsizes(j), "Index out of bounds");
74                                        v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j);
75                                }
76                                return v;
77                        }
78
79                        //! convert dimension indeces into linear index, the indexing is in the same way as in \c next_vec()
80                        long linear_index (const ivec inds){
81                                long ind=0;
82                                bdm_assert_debug(inds.length() == dim, "Improper indices inds");
83                                int dim_skips=1; // skips in active dimension, in the first dimension, the skips are 1 index per value
84                                               
85                                for (int j=0; j<dim; j++){ 
86                                        ind += dim_skips*(inds(j)); // add shift in linear index caused by this dimension
87                                        dim_skips*=gridsizes(j);  // indeces in the next dimension are repeated with period gridsizes(j) times greater that in this dimesion
88                                }
89                                return ind;
90                        } 
91                        //! set the first vector to corner and store result in actvec
92                        const vec& first_vec(){
93                                for ( int j = 0; j < dim; j++ ) {
94                                        actvec ( j ) = ranges(j)(0) + 0.5*steps(j);
95                                        actvec_ind(j) = 0;
96                                }
97                                return actvec;
98                        }
99                        //! Get next active vector, call ONLY after first_vector()!
100                        const vec& next_vec() {
101                                // go through all dimensions
102                                int j=0;
103                                while (j<dim) {
104                                        if ( actvec_ind ( j ) == gridsizes ( j ) - 1 ) { //j-th index is full
105                                                actvec_ind ( j ) = 0; //shift back
106                                                actvec ( j ) = ranges ( j ) ( 0 ) + 0.5*steps(j);
107                                                j++;
108                                        } else {
109                                                actvec_ind ( j ) ++;
110                                                actvec ( j ) += steps ( j );
111                                                break;
112                                        }
113                                }
114                                return actvec;
115                        } 
116                       
117                        ivec nearest_point(const vec &val){
118                                ivec inds;
119                                inds.set_size(dim);
120                                for(int j=0; j<dim; j++){
121                                        if (val(j)<ranges(j)(0))
122                                                inds(j) = 0;
123                                        else {
124                                                if (val(j)>ranges(j)(1))
125                                                        inds(j) = gridsizes(j)-1;
126                                                else {
127                                                        inds(j) = ::round(val(j) - ranges(j)(0)/ steps(j));
128                                                }
129                                        }
130                                }
131                                return inds;
132                        }
133
134                        //! Access function
135                        int points() const {return Npoints;}
136                       
137                        //! access function
138                        const vec& _steps() const {return steps;}
139                       
140                        void from_setting (const Setting &set) {
141                                UI::get (ranges , set, "ranges", UI::compulsory);
142                                UI::get (gridsizes, set, "gridsizes", UI::compulsory);
143                                initialize();
144                        }
145        };
146        UIREGISTER(rectangular_support);
147       
148        //! Discrete support with stored support points
149        class discrete_support: public root{
150                protected:
151                        //! storage of support points
152                        Array<vec> Spoints;
153                        //! index in iterators
154                        int idx;
155                public:
156                        //! Default constructor
157                        discrete_support() : Spoints(0), idx(0){}
158                        //! Access function
159                        int points() const {return Spoints.length();}
160                        //! set the first vector to corner and store result in actvec
161                        const vec& first_vec(){bdm_assert_debug(Spoints.length()>0,"Empty support");idx=0; return Spoints(idx);}
162                        //! set next vector after calling first_vec()
163                        const vec& next_vec(){bdm_assert_debug(Spoints.length()>idx-1,"Out of support points"); return Spoints(++idx);}
164                       
165                        /*!
166                        \code
167                          class = "discrete_support";
168                          points = ( [1,2..], [2,2..], ...); // list of points
169                           === OR ===
170                          epdf = {class="epdf_offspring",...}; // epdf rfom which to sample
171                          npoints = 100;                     // number of samples
172                        \endcode
173                        */
174                        void from_setting (const Setting &set) {
175                                UI::get (Spoints, set, "points", UI::compulsory);
176                                if (points()<1){
177                                        int npoints;
178                                        shared_ptr<epdf> ep= UI::build<epdf>(set, "epdf", UI::compulsory);
179                                        if (!UI::get(npoints,set,"npoints",UI::optional)){npoints=100;}
180                                       
181                                        //sample
182                                        Spoints.set_size(npoints);
183                                        for(int i=0; i<points(); i++){Spoints(i)=ep->sample();}
184                                }
185                        }
186                        //! access function
187                        Array<vec> & _Spoints() {return Spoints;}
188        };
189        UIREGISTER(discrete_support);
190       
191        class grid_fnc: public fnc{
192                protected:
193                        rectangular_support sup;
194                        vec values;
195                public:
196                        //! constructor function
197                        void set_support(rectangular_support &sup0){sup=sup0; values=zeros(sup.points());}
198                        //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
199                        void set_values(double (*evalptr)(const vec&)){
200                                if (sup.points()>0){
201                                        values(0) = (*evalptr)(sup.first_vec());
202                                        for (int j=1; j<sup.points(); j++){ 
203                                                values(j)=(*evalptr)(sup.next_vec());
204                                        }
205                                }
206                        } 
207                        //! get value nearest to the given point
208                        double nearest_val(const vec &val){return values( sup.linear_index(sup.nearest_point(val)));}
209                       
210                        vec eval(const vec &val){return vec_1(nearest_val(val));}
211        };
212        UIREGISTER(grid_fnc);
213
214        //! Piecewise constant pdf on rectangular support
215        //! Each point on the grid represents a centroid around which the density is constant.
216        //! This is a trivial point-mass density where all points have the same mass.
217        class egrid: public epdf{
218                protected:
219                        rectangular_support sup;
220                        vec values;
221                public:
222                        //! we assume that evallog is not called too often otherwise we should cache log(values)
223                        double evallog(const vec &val){return log(values( sup.linear_index(sup.nearest_point(val))));}
224                       
225        };
226}
227#endif //DISCR_H
Note: See TracBrowser for help on using the browser.