| 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(gridsizes.length() == dim, "Incompatible dimensions of input"); |
| 58 | | Npoints = prod(gridsizes); |
| 59 | | bdm_assert(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 | | } |
| | 23 | //! Rectangular support |
| | 24 | //! Support points are located inbetween ranges! For example: |
| | 25 | //! For ranges=[0,1] and gridsizes=[1] the support point is 0.5 |
| | 26 | class rectangular_support: public root { |
| | 27 | protected: |
| | 28 | //! Array of boundaries (2D vectors: [begining,end]) for each dimension |
| | 29 | Array<vec> ranges; |
| | 30 | //! Number of support points in each dimension |
| | 31 | ivec gridsizes; |
| | 32 | //! dimension |
| | 33 | int dim; |
| | 34 | //! Number of data points |
| | 35 | int Npoints; |
| | 36 | //! active vector for first_vec and next_vec |
| | 37 | vec actvec; |
| | 38 | //! indeces of active vector |
| | 39 | vec actvec_ind; |
| | 40 | //! length of steps in each dimension |
| | 41 | vec steps; |
| | 42 | public: |
| | 43 | //! default constructor |
| | 44 | rectangular_support() : dim ( 0 ), Npoints ( 0 ) { |
| | 45 | } |
| 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 | | } |
| | 47 | //! set parameters |
| | 48 | void set_parameters ( const Array<vec> &ranges0, const ivec &gridsize0 ); |
| 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); |
| | 50 | //! Internal functio to set temporaries correctly |
| | 51 | void initialize(); |
| 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 | | }; |
| | 53 | //! return vector at position given by vector of indeces |
| | 54 | vec get_vec ( const ivec &inds ); |
| | 55 | |
| | 56 | //! convert dimension indeces into linear index, the indexing is in the same way as in \c next_vec() |
| | 57 | long linear_index ( const ivec inds ); |
| | 58 | |
| | 59 | //! set the first vector to corner and store result in actvec |
| | 60 | const vec& first_vec(); |
| | 61 | |
| | 62 | //! Get next active vector, call ONLY after first_vector()! |
| | 63 | const vec& next_vec(); |
| | 64 | |
| | 65 | //! \todo to je asi navic .. v predkovi! |
| | 66 | ivec nearest_point ( const vec &val ); |
| | 67 | |
| | 68 | //! Access function |
| | 69 | int points() const { |
| | 70 | return Npoints; |
| | 71 | } |
| | 72 | |
| | 73 | //! access function |
| | 74 | //! \todo opet pouze do potomka.. |
| | 75 | const vec& _steps() const { |
| | 76 | return steps; |
| | 77 | } |
| | 78 | |
| | 79 | void from_setting ( const Setting &set ); |
| | 80 | }; |
| | 81 | UIREGISTER ( rectangular_support ); |
| | 82 | |
| | 83 | //! Discrete support with stored support points |
| | 84 | class discrete_support: public root { |
| | 85 | protected: |
| | 86 | //! storage of support points |
| | 87 | Array<vec> Spoints; |
| | 88 | //! index in iterators |
| | 89 | int idx; |
| | 90 | public: |
| | 91 | //! Default constructor |
| | 92 | discrete_support() : Spoints ( 0 ), idx ( 0 ) {} |
| | 93 | //! Access function |
| | 94 | int points() const { |
| | 95 | return Spoints.length(); |
| | 96 | } |
| | 97 | //! set the first vector to corner and store result in actvec |
| | 98 | const vec& first_vec() { |
| | 99 | bdm_assert_debug ( Spoints.length() > 0, "Empty support" ); |
| | 100 | idx = 0; |
| | 101 | return Spoints ( idx ); |
| | 102 | } |
| | 103 | //! set next vector after calling first_vec() |
| | 104 | const vec& next_vec() { |
| | 105 | bdm_assert_debug ( Spoints.length() > idx - 1, "Out of support points" ); |
| | 106 | return Spoints ( ++idx ); |
| | 107 | } |
| | 108 | |
| | 109 | /*! |
| | 110 | \code |
| | 111 | class = "discrete_support"; |
| | 112 | points = ( [1,2..], [2,2..], ...); // list of points |
| | 113 | === OR === |
| | 114 | epdf = {class="epdf_offspring",...}; // epdf rfom which to sample |
| | 115 | npoints = 100; // number of samples |
| | 116 | \endcode |
| | 117 | */ |
| | 118 | void from_setting ( const Setting &set ); |
| | 119 | |
| | 120 | //! access function |
| | 121 | Array<vec> & _Spoints() { |
| | 122 | return Spoints; |
| | 123 | } |
| | 124 | }; |
| | 125 | UIREGISTER ( discrete_support ); |
| | 126 | |
| | 127 | class grid_fnc: public fnc { |
| | 128 | protected: |
| | 129 | rectangular_support sup; |
| | 130 | vec values; |
| | 131 | public: |
| | 132 | //! constructor function |
| | 133 | void set_support ( rectangular_support &sup0 ) { |
| | 134 | sup = sup0; |
| | 135 | values = zeros ( sup.points() ); |
| | 136 | } |
| | 137 | //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer |
| | 138 | void set_values ( double ( *evalptr ) ( const vec& ) ); |
| | 139 | |
| | 140 | //! get value nearest to the given point |
| | 141 | double nearest_val ( const vec &val ) { |
| | 142 | return values ( sup.linear_index ( sup.nearest_point ( val ) ) ); |
| | 143 | } |
| | 144 | |
| | 145 | vec eval ( const vec &val ) { |
| | 146 | return vec_1 ( nearest_val ( val ) ); |
| | 147 | } |
| | 148 | }; |
| | 149 | UIREGISTER ( grid_fnc ); |
| | 150 | |
| | 151 | //! Piecewise constant pdf on rectangular support |
| | 152 | //! Each point on the grid represents a centroid around which the density is constant. |
| | 153 | //! This is a trivial point-mass density where all points have the same mass. |
| | 154 | class egrid: public epdf { |
| | 155 | protected: |
| | 156 | rectangular_support sup; |
| | 157 | vec values; |
| | 158 | public: |
| | 159 | //! we assume that evallog is not called too often otherwise we should cache log(values) |
| | 160 | double evallog ( const vec &val ) { |
| | 161 | return log ( values ( sup.linear_index ( sup.nearest_point ( val ) ) ) ); |
| | 162 | } |
| | 163 | |
| | 164 | }; |
