[556] | 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 | |
---|
| 21 | namespace bdm |
---|
| 22 | { |
---|
| 23 | |
---|
| 24 | //! Rectangular support |
---|
| 25 | //! Support ponits 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(){}; |
---|
| 46 | |
---|
| 47 | //! set parameters |
---|
| 48 | void set_parameters(const Array<vec> &ranges0, const ivec &gridsize0){ |
---|
| 49 | ranges=ranges0; |
---|
| 50 | gridsizes=gridsize0; |
---|
| 51 | initialize(); |
---|
| 52 | } |
---|
| 53 | //! Internal functio to set temporaries correctly |
---|
| 54 | void initialize() { |
---|
| 55 | dim = ranges.length(); |
---|
| 56 | it_assert_debug(gridsizes.length()==dim,"Incompatible dimensions of input"); |
---|
| 57 | Npoints = prod(gridsizes); |
---|
| 58 | it_assert_debug(Npoints>0,"Wrong input parameters"); |
---|
| 59 | |
---|
| 60 | //precompute steps |
---|
[557] | 61 | steps.set_size(dim); |
---|
[556] | 62 | for ( int j = 0; j < dim; j++ ) { |
---|
| 63 | steps ( j ) = ( ranges ( j ) ( 1 ) - ranges(j)(0) ) / gridsizes ( j ); |
---|
| 64 | } |
---|
| 65 | actvec.set_size(dim); |
---|
| 66 | actvec_ind.set_size(dim); |
---|
| 67 | } |
---|
[558] | 68 | //! return vector at position given by vector of indeces |
---|
| 69 | vec get_vec(const ivec &inds){ |
---|
| 70 | vec v; |
---|
| 71 | for ( int j = 0; j < dim; j++ ) { |
---|
| 72 | it_assert_debug(inds(j)<gridsizes(j), "Index out of bounds"); |
---|
| 73 | v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j); |
---|
| 74 | } |
---|
| 75 | return v; |
---|
| 76 | } |
---|
[556] | 77 | |
---|
[558] | 78 | //! convert dimension indeces into linear index, the indexing is in the same way as in \c next_vec() |
---|
| 79 | long linear_index (const ivec inds){ |
---|
| 80 | long ind=0; |
---|
| 81 | it_assert_debug(inds.length()==dim,"Improper indeces inds"); |
---|
| 82 | int dim_skips=1; // skips in active dimension, in the first dimension, the skips are 1 index per value |
---|
| 83 | |
---|
| 84 | for (int j=0; j<dim; j++){ |
---|
| 85 | ind += dim_skips*(inds(j)); // add shift in linear index caused by this dimension |
---|
| 86 | dim_skips*=gridsizes(j); // indeces in the next dimension are repeated with period gridsizes(j) times greater that in this dimesion |
---|
| 87 | } |
---|
| 88 | return ind; |
---|
| 89 | } |
---|
[556] | 90 | //! set the first corner to actvec |
---|
| 91 | const vec& first_vec(){ |
---|
| 92 | for ( int j = 0; j < dim; j++ ) { |
---|
| 93 | actvec ( j ) = ranges(j)(0) + 0.5*steps(j); |
---|
| 94 | actvec_ind(j) = 0; |
---|
| 95 | } |
---|
| 96 | return actvec; |
---|
| 97 | } |
---|
| 98 | //! Get next active vector, call ONLY after first_vector()! |
---|
| 99 | const vec& next_vec() { |
---|
| 100 | // go through all dimensions |
---|
[558] | 101 | int j=0; |
---|
| 102 | while (j<dim) { |
---|
[556] | 103 | if ( actvec_ind ( j ) == gridsizes ( j ) - 1 ) { //j-th index is full |
---|
| 104 | actvec_ind ( j ) = 0; //shift back |
---|
[557] | 105 | actvec ( j ) = ranges ( j ) ( 0 ) + 0.5*steps(j); |
---|
[558] | 106 | j++; |
---|
[556] | 107 | } else { |
---|
| 108 | actvec_ind ( j ) ++; |
---|
| 109 | actvec ( j ) += steps ( j ); |
---|
| 110 | break; |
---|
| 111 | } |
---|
| 112 | } |
---|
| 113 | return actvec; |
---|
[557] | 114 | } |
---|
[558] | 115 | |
---|
| 116 | ivec nearest_point(const vec &val){ |
---|
| 117 | ivec inds; |
---|
| 118 | inds.set_size(dim); |
---|
| 119 | for(int j=0; j<dim; j++){ |
---|
| 120 | if (val(j)<ranges(j)(0)) |
---|
| 121 | inds(j) = 0; |
---|
| 122 | else { |
---|
| 123 | if (val(j)>ranges(j)(1)) |
---|
| 124 | inds(j) = gridsizes(j)-1; |
---|
| 125 | else { |
---|
| 126 | inds(j) = ::round(val(j) - ranges(j)(0)/ steps(j)); |
---|
| 127 | } |
---|
| 128 | } |
---|
| 129 | } |
---|
| 130 | return inds; |
---|
| 131 | } |
---|
| 132 | |
---|
[557] | 133 | //! Access function |
---|
[556] | 134 | int points() const {return Npoints;} |
---|
| 135 | |
---|
[557] | 136 | //! access function |
---|
| 137 | const vec& _steps() const {return steps;} |
---|
| 138 | |
---|
[556] | 139 | void from_setting (const Setting &set) { |
---|
| 140 | UI::get (ranges , set, "ranges", UI::compulsory); |
---|
| 141 | UI::get (gridsizes, set, "gridsizes", UI::compulsory); |
---|
| 142 | initialize(); |
---|
| 143 | } |
---|
| 144 | }; |
---|
| 145 | UIREGISTER(rectangular_support); |
---|
[558] | 146 | |
---|
| 147 | class grid_fnc: public fnc{ |
---|
| 148 | protected: |
---|
| 149 | rectangular_support sup; |
---|
| 150 | vec values; |
---|
| 151 | public: |
---|
| 152 | //! constructor function |
---|
| 153 | void set_support(rectangular_support &sup0){sup=sup0; values=zeros(sup.points());} |
---|
| 154 | //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer |
---|
| 155 | void set_values(double (*evalptr)(const vec&)){ |
---|
| 156 | if (sup.points()>0){ |
---|
| 157 | values(0) = (*evalptr)(sup.first_vec()); |
---|
| 158 | for (int j=1; j<sup.points(); j++){ |
---|
| 159 | values(j)=(*evalptr)(sup.next_vec()); |
---|
| 160 | } |
---|
| 161 | } |
---|
| 162 | } |
---|
| 163 | //! get value nearest to the given point |
---|
| 164 | double nearest_val(const vec &val){return values( sup.linear_index(sup.nearest_point(val)));} |
---|
| 165 | |
---|
| 166 | vec eval(const vec &val){return vec_1(nearest_val(val));} |
---|
| 167 | }; |
---|
| 168 | UIREGISTER(grid_fnc); |
---|
| 169 | |
---|
| 170 | //! Piecewise constant pdf on rectangular support |
---|
| 171 | //! Each point on the grid represents a centroid around which the density is constant. |
---|
| 172 | //! This is a trivial point-mass density where all points have the same mass. |
---|
| 173 | class egrid: public epdf{ |
---|
| 174 | protected: |
---|
| 175 | rectangular_support sup; |
---|
| 176 | vec values; |
---|
| 177 | public: |
---|
| 178 | //! we assume that evallog is not called too often otherwise we should cache log(values) |
---|
| 179 | double evallog(const vec &val){return log(values( sup.linear_index(sup.nearest_point(val))));} |
---|
| 180 | |
---|
| 181 | }; |
---|
[556] | 182 | } |
---|
| 183 | #endif //DISCR_H |
---|