Changeset 558

Timestamp:

08/19/09 01:41:25 (15 years ago)

Author:

smidl

Message:

New stuff for discrete grid world - pdf and function

Files:

• library/bdm/stat/discrete.h (modified) (4 diffs)

library/bdm/stat/discrete.h

@@ -66 +66 @@
                                 actvec_ind.set_size(dim);
                         }
+                        //! return vector at position given by vector of indeces
+                        vec get_vec(const ivec &inds){
+                                vec v;
+                                for ( int j = 0; j < dim; j++ ) {
+                                        it_assert_debug(inds(j)<gridsizes(j), "Index out of bounds");
+                                        v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j);
+                                }
+                                return v;
+                        }
 
+                        //! convert dimension indeces into linear index, the indexing is in the same way as in \c next_vec()
+                        long linear_index (const ivec inds){
+                                long ind=0;
+                                it_assert_debug(inds.length()==dim,"Improper indeces inds");
+                                int dim_skips=1; // skips in active dimension, in the first dimension, the skips are 1 index per value
+                                                
+                                for (int j=0; j<dim; j++){ 
+                                        ind += dim_skips*(inds(j)); // add shift in linear index caused by this dimension
+                                        dim_skips*=gridsizes(j);  // indeces in the next dimension are repeated with period gridsizes(j) times greater that in this dimesion
+                                }
+                                return ind;
+                        } 
                         //! set the first corner to actvec
                         const vec& first_vec(){
@@ -78 +99 @@
                         const vec& next_vec() {
                                 // go through all dimensions
-                                for ( int j = 0; j < dim; j++ ) {
+                                int j=0;
+                                while (j<dim) {
                                         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;
-                                                        }
-
+                                                j++;
                                         } else {
                                                 actvec_ind ( j ) ++;
@@ -99 +113 @@
                                 return actvec;
                         } 
+                        
+                        ivec nearest_point(const vec &val){
+                                ivec inds;
+                                inds.set_size(dim);
+                                for(int j=0; j<dim; j++){
+                                        if (val(j)<ranges(j)(0))
+                                                inds(j) = 0;
+                                        else {
+                                                if (val(j)>ranges(j)(1))
+                                                        inds(j) = gridsizes(j)-1;
+                                                else {
+                                                        inds(j) = ::round(val(j) - ranges(j)(0)/ steps(j));
+                                                }
+                                        }
+                                }
+                                return inds;
+                        }
+
                         //! Access function
                         int points() const {return Npoints;}
@@ -112 +144 @@
         };
         UIREGISTER(rectangular_support);
+        
+        class grid_fnc: public fnc{
+                protected:
+                        rectangular_support sup;
+                        vec values;
+                public:
+                        //! constructor function
+                        void set_support(rectangular_support &sup0){sup=sup0; values=zeros(sup.points());}
+                        //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
+                        void set_values(double (*evalptr)(const vec&)){
+                                if (sup.points()>0){
+                                        values(0) = (*evalptr)(sup.first_vec());
+                                        for (int j=1; j<sup.points(); j++){ 
+                                                values(j)=(*evalptr)(sup.next_vec());
+                                        }
+                                }
+                        } 
+                        //! get value nearest to the given point
+                        double nearest_val(const vec &val){return values( sup.linear_index(sup.nearest_point(val)));}
+                        
+                        vec eval(const vec &val){return vec_1(nearest_val(val));}
+        };
+        UIREGISTER(grid_fnc);
+
+        //! Piecewise constant pdf on rectangular support
+        //! Each point on the grid represents a centroid around which the density is constant.
+        //! This is a trivial point-mass density where all points have the same mass.
+        class egrid: public epdf{
+                protected:
+                        rectangular_support sup;
+                        vec values;
+                public:
+                        //! we assume that evallog is not called too often otherwise we should cache log(values)
+                        double evallog(const vec &val){return log(values( sup.linear_index(sup.nearest_point(val))));}
+                        
+        };
 }
 #endif //DISCR_H