- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# discrete.h

Go to the documentation of this file.

    00001
    00013 #ifndef DISCR_H
    00014 #define DISCR_H
    00015
    00016
    00017 #include "../shared_ptr.h"
    00018 #include "../base/bdmbase.h"
    00019 #include "../math/chmat.h"
    00020
    00021 namespace bdm
    00022 {
    00023
    00027         class rectangular_support: public root {
    00028                 protected:
    00030                         Array<vec> ranges;
    00032                         ivec gridsizes;
    00034                         int dim;
    00036                         int Npoints;
    00038                         vec actvec;
    00040                         vec actvec_ind;
    00042                         vec steps;
    00043                 public:
    00045                         rectangular_support() : dim(0), Npoints(0) {
    00046                         }
    00047
    00049                         void set_parameters(const Array<vec> &ranges0, const ivec &gridsize0){
    00050                                 ranges=ranges0;
    00051                                 gridsizes=gridsize0;
    00052                                 initialize();
    00053                         }
    00055                         void initialize() {
    00056                                 dim = ranges.length();
    00057                                 bdm_assert_debug(gridsizes.length() == dim, "Incompatible dimensions of input");
    00058                                 Npoints = prod(gridsizes);
    00059                                 bdm_assert_debug(Npoints > 0, "Wrong input parameters");
    00060
    00061                                 //precompute steps
    00062                                 steps.set_size(dim);
    00063                                 for ( int j = 0; j < dim; j++ ) {
    00064                                         steps ( j ) = ( ranges ( j ) ( 1 ) - ranges(j)(0) ) / gridsizes ( j );
    00065                                 }
    00066                                 actvec.set_size(dim);
    00067                                 actvec_ind.set_size(dim);
    00068                         }
    00070                         vec get_vec(const ivec &inds){
    00071                                 vec v ( dim );
    00072                                 for ( int j = 0; j < dim; j++ ) {
    00073                                         bdm_assert_debug(inds(j) < gridsizes(j), "Index out of bounds");
    00074                                         v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j);
    00075                                 }
    00076                                 return v;
    00077                         }
    00078
    00080                         long linear_index (const ivec inds){
    00081                                 long ind=0;
    00082                                 bdm_assert_debug(inds.length() == dim, "Improper indices inds");
    00083                                 int dim_skips=1; // skips in active dimension, in the first dimension, the skips are 1 index per value
    00084
    00085                                 for (int j=0; j<dim; j++){
    00086                                         ind += dim_skips*(inds(j)); // add shift in linear index caused by this dimension
    00087                                         dim_skips*=gridsizes(j);  // indeces in the next dimension are repeated with period gridsizes(j) times greater that in this dimesion
    00088                                 }
    00089                                 return ind;
    00090                         }
    00092                         const vec& first_vec(){
    00093                                 for ( int j = 0; j < dim; j++ ) {
    00094                                         actvec ( j ) = ranges(j)(0) + 0.5*steps(j);
    00095                                         actvec_ind(j) = 0;
    00096                                 }
    00097                                 return actvec;
    00098                         }
    00100                         const vec& next_vec() {
    00101                                 // go through all dimensions
    00102                                 int j=0;
    00103                                 while (j<dim) {
    00104                                         if ( actvec_ind ( j ) == gridsizes ( j ) - 1 ) { //j-th index is full
    00105                                                 actvec_ind ( j ) = 0; //shift back
    00106                                                 actvec ( j ) = ranges ( j ) ( 0 ) + 0.5*steps(j);
    00107                                                 j++;
    00108                                         } else {
    00109                                                 actvec_ind ( j ) ++;
    00110                                                 actvec ( j ) += steps ( j );
    00111                                                 break;
    00112                                         }
    00113                                 }
    00114                                 return actvec;
    00115                         }
    00116
    00117                         ivec nearest_point(const vec &val){
    00118                                 ivec inds;
    00119                                 inds.set_size(dim);
    00120                                 for(int j=0; j<dim; j++){
    00121                                         if (val(j)<ranges(j)(0))
    00122                                                 inds(j) = 0;
    00123                                         else {
    00124                                                 if (val(j)>ranges(j)(1))
    00125                                                         inds(j) = gridsizes(j)-1;
    00126                                                 else {
    00127                                                         inds(j) = ::round(val(j) - ranges(j)(0)/ steps(j));
    00128                                                 }
    00129                                         }
    00130                                 }
    00131                                 return inds;
    00132                         }
    00133
    00135                         int points() const {return Npoints;}
    00136
    00138                         const vec& _steps() const {return steps;}
    00139
    00140                         void from_setting (const Setting &set) {
    00141                                 UI::get (ranges , set, "ranges", UI::compulsory);
    00142                                 UI::get (gridsizes, set, "gridsizes", UI::compulsory);
    00143                                 initialize();
    00144                         }
    00145         };
    00146         UIREGISTER(rectangular_support);
    00147
    00149         class discrete_support: public root{
    00150                 protected:
    00152                         Array<vec> Spoints;
    00154                         int idx;
    00155                 public:
    00157                         discrete_support() : Spoints(0), idx(0){}
    00159                         int points() const {return Spoints.length();}
    00161                         const vec& first_vec(){bdm_assert_debug(Spoints.length()>0,"Empty support");idx=0; return Spoints(idx);}
    00163                         const vec& next_vec(){bdm_assert_debug(Spoints.length()>idx-1,"Out of support points"); return Spoints(++idx);}
    00164
    00174                         void from_setting (const Setting &set) {
    00175                                 UI::get (Spoints, set, "points", UI::compulsory);
    00176                                 if (points()<1){
    00177                                         int npoints;
    00178                                         shared_ptr<epdf> ep= UI::build<epdf>(set, "epdf", UI::compulsory);
    00179                                         if (!UI::get(npoints,set,"npoints",UI::optional)){npoints=100;}
    00180
    00181                                         //sample
    00182                                         Spoints.set_size(npoints);
    00183                                         for(int i=0; i<points(); i++){Spoints(i)=ep->sample();}
    00184                                 }
    00185                         }
    00187                         Array<vec> & _Spoints() {return Spoints;}
    00188         };
    00189         UIREGISTER(discrete_support);
    00190
    00191         class grid_fnc: public fnc{
    00192                 protected:
    00193                         rectangular_support sup;
    00194                         vec values;
    00195                 public:
    00197                         void set_support(rectangular_support &sup0){sup=sup0; values=zeros(sup.points());}
    00199                         void set_values(double (*evalptr)(const vec&)){
    00200                                 if (sup.points()>0){
    00201                                         values(0) = (*evalptr)(sup.first_vec());
    00202                                         for (int j=1; j<sup.points(); j++){
    00203                                                 values(j)=(*evalptr)(sup.next_vec());
    00204                                         }
    00205                                 }
    00206                         }
    00208                         double nearest_val(const vec &val){return values( sup.linear_index(sup.nearest_point(val)));}
    00209
    00210                         vec eval(const vec &val){return vec_1(nearest_val(val));}
    00211         };
    00212         UIREGISTER(grid_fnc);
    00213
    00217         class egrid: public epdf{
    00218                 protected:
    00219                         rectangular_support sup;
    00220                         vec values;
    00221                 public:
    00223                         double evallog(const vec &val){return log(values( sup.linear_index(sup.nearest_point(val))));}
    00224
    00225         };
    00226 }
    00227 #endif //DISCR_H