discrete.h

Go to the documentation of this file.

#ifndef DISCR_H
#define DISCR_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm
{

        class rectangular_support: public root {
                protected:
                        Array<vec> ranges;
                        ivec gridsizes;
                        int dim;
                        int Npoints;
                        vec actvec;
                        vec actvec_ind;
                        vec steps;
                public:
                        rectangular_support() : dim(0), Npoints(0) {
                        }
                        
                        void set_parameters(const Array<vec> &ranges0, const ivec &gridsize0){
                                ranges=ranges0;
                                gridsizes=gridsize0;
                                initialize();
                        }
                        void initialize() {
                                dim = ranges.length();
                                bdm_assert(gridsizes.length() == dim, "Incompatible dimensions of input");
                                Npoints = prod(gridsizes);
                                bdm_assert(Npoints > 0, "Wrong input parameters");
                                
                                //precompute steps
                                steps.set_size(dim);
                                for ( int j = 0; j < dim; j++ ) {
                                        steps ( j ) = ( ranges ( j ) ( 1 ) - ranges(j)(0) ) / gridsizes ( j );
                                }
                                actvec.set_size(dim);
                                actvec_ind.set_size(dim);
                        }
                        vec get_vec(const ivec &inds){
                                vec v ( dim );
                                for ( int j = 0; j < dim; j++ ) {
                                        bdm_assert_debug(inds(j) < gridsizes(j), "Index out of bounds");
                                        v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j);
                                }
                                return v;
                        }

                        long linear_index (const ivec inds){
                                long ind=0;
                                bdm_assert_debug(inds.length() == dim, "Improper indices 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;
                        } 
                        const vec& first_vec(){
                                for ( int j = 0; j < dim; j++ ) {
                                        actvec ( j ) = ranges(j)(0) + 0.5*steps(j);
                                        actvec_ind(j) = 0;
                                }
                                return actvec;
                        }
                        const vec& next_vec() {
                                // go through all dimensions
                                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);
                                                j++;
                                        } else {
                                                actvec_ind ( j ) ++;
                                                actvec ( j ) += steps ( j );
                                                break;
                                        }
                                }
                                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;
                        }

                        int points() const {return Npoints;}
                        
                        const vec& _steps() const {return steps;}
                        
                        void from_setting (const Setting &set) {
                                UI::get (ranges , set, "ranges", UI::compulsory);
                                UI::get (gridsizes, set, "gridsizes", UI::compulsory);
                                initialize();
                        }
        };
        UIREGISTER(rectangular_support);
        
        class discrete_support: public root{
                protected:
                        Array<vec> Spoints;
                        int idx;
                public:
                        discrete_support() : Spoints(0), idx(0){}
                        int points() const {return Spoints.length();}
                        const vec& first_vec(){bdm_assert_debug(Spoints.length()>0,"Empty support");idx=0; return Spoints(idx);}
                        const vec& next_vec(){bdm_assert_debug(Spoints.length()>idx-1,"Out of support points"); return Spoints(++idx);}
                        
                        void from_setting (const Setting &set) {
                                UI::get (Spoints, set, "points", UI::compulsory);
                                if (points()<1){
                                        int npoints;
                                        shared_ptr<epdf> ep= UI::build<epdf>(set, "epdf", UI::compulsory);
                                        if (!UI::get(npoints,set,"npoints",UI::optional)){npoints=100;}
                                        
                                        //sample
                                        Spoints.set_size(npoints);
                                        for(int i=0; i<points(); i++){Spoints(i)=ep->sample();}
                                }
                        }
                        Array<vec> & _Spoints() {return Spoints;}
        };
        UIREGISTER(discrete_support);
        
        class grid_fnc: public fnc{
                protected:
                        rectangular_support sup;
                        vec values;
                public:
                        void set_support(rectangular_support &sup0){sup=sup0; values=zeros(sup.points());}
                        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());
                                        }
                                }
                        } 
                        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);

        class egrid: public epdf{
                protected:
                        rectangular_support sup;
                        vec values;
                public:
                        double evallog(const vec &val){return log(values( sup.linear_index(sup.nearest_point(val))));}
                        
        };
}
#endif //DISCR_H

---