Changeset 1188

Timestamp:

09/16/10 13:21:59 (14 years ago)

Author:

smidl

Message:

support has abstract base

Location:

library/bdm/stat

Files:

• discrete.cpp (modified) (4 diffs)
• discrete.h (modified) (9 diffs)

library/bdm/stat/discrete.cpp

@@ -50 +50 @@
         actvec_ind ( j ) = 0;
     }
+    actvolume = prod(steps);
     return actvec;
 }
@@ -66 +67 @@
         }
     }
+    // no need to set active volume
     return actvec;
 }
@@ -111 +113 @@
         }
     }
+    
+    UI::get (volumes, set, "volumes",UI::optional);
+        if (volumes.length()<1){
+                volumes = 1.0/points()*ones(points());
+        } 
+        
 }
 
 void grid_fnc::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() );
+    if ( sup->points() > 0 ) {
+        values ( 0 ) = ( *evalptr ) ( sup->first_vec() );
+        for ( int j = 1; j < sup->points(); j++ ) {
+            values ( j ) = ( *evalptr ) ( sup->next_vec() );
         }
     }
@@ -123 +131 @@
 
 void grid_fnc::set_values ( const epdf &ep ) {
-    if ( sup.points() > 0 ) {
-        values ( 0 ) =  exp ( ep.evallog ( sup.first_vec() ) );
-        for ( int j = 1; j < sup.points(); j++ ) {
-            values ( j ) = exp ( ep.evallog ( sup.next_vec() ) );
+    if ( sup->points() > 0 ) {
+        values ( 0 ) =  exp ( ep.evallog ( sup->first_vec() ) );
+        for ( int j = 1; j < sup->points(); j++ ) {
+            values ( j ) = exp ( ep.evallog ( sup->next_vec() ) );
         }
     }
 }
 
+double grid_fnc::integrate ( ) {
+        double intgr=0.0;
+        if ( sup->points() > 0 ) {
+                sup->first_vec(); // also first volume
+                intgr+= values(0) * sup->act_volume();
+                for ( int j = 1; j < sup->points(); j++ ) {
+                        sup->next_vec();
+                        intgr+= values(j)*sup->act_volume();
+                }
+        }
+        return intgr;
 }
+
+}

library/bdm/stat/discrete.h

@@ -21 +21 @@
 namespace bdm {
 
+        class support_base: public root {
+        protected:
+                //! dimension
+                int dim;
+                //! Number of data points
+                int Npoints;
+                //! active vector for first_vec and next_vec
+                vec actvec;
+                //! volume of the active point
+                double actvolume;
+        public:
+                //! default constructor
+                support_base() : dim ( 0 ), Npoints ( 0 ) {
+                }
+                                                
+                //! set the first vector to corner and store result in actvec
+                virtual const vec& first_vec() NOT_IMPLEMENTED(empty_vec);
+                
+                //! Get next active vector, call ONLY after first_vector()!
+                virtual const vec& next_vec() NOT_IMPLEMENTED(empty_vec);;
+                
+                //! return active vector, call ONLY after first_vector()!
+                virtual const vec& act_vec() {
+                        return actvec;
+                };
+                //! return active vector, call ONLY after first_vector()!
+                virtual const double& act_volume() {
+                        return actvolume;
+                };
+                
+                //! Access function
+                int points() const {
+                        return Npoints;
+                }
+                                
+        };
+        UIREGISTER ( support_base );
+        
 //! Rectangular support
 //! Support points are located inbetween ranges! For example:
 //! For ranges=[0,1] and gridsizes=[1] the support point is 0.5
-class rectangular_support: public root {
+class rectangular_support: public support_base {
 protected:
     //! Array of boundaries (2D vectors: [begining,end]) for each dimension
@@ -30 +68 @@
     //! Number of support points in each dimension
     ivec gridsizes;
-    //! dimension
-    int dim;
-    //! Number of data points
-    int Npoints;
-    //! active vector for first_vec and next_vec
-    vec actvec;
-    //! indices of active vector
-    vec actvec_ind;
-    //! length of steps in each dimension
-    vec steps;
+        //! indices of active vector
+        vec actvec_ind;
+        //
+        vec steps;
 public:
     //! default constructor
-    rectangular_support() : dim ( 0 ), Npoints ( 0 ) {
+    rectangular_support() : support_base() {
     }
 
@@ -87 +119 @@
 
 //! Discrete support with stored support points
-class discrete_support: public root {
+class discrete_support: public support_base {
 protected:
     //! storage of support points
@@ -93 +125 @@
     //! index in iterators
     int idx;
+        //! volumes for each point
+        vec volumes;
 public:
     //! Default constructor
-    discrete_support() : Spoints ( 0 ), idx ( 0 ) {}
+    discrete_support() : Spoints ( 0 ), idx ( 0 ), volumes(0) {}
     //! Access function
     int points() const {
@@ -102 +136 @@
     //! set the first vector to corner and store result in actvec
     const vec& first_vec() {
-        bdm_assert_debug ( Spoints.length() > 0, "Empty support" );
+
+                bdm_assert_debug ( Spoints.length() > 0, "Empty support" );
         idx = 0;
+                actvolume=volumes(0);
         return Spoints ( idx );
     }
@@ -109 +145 @@
     const vec& next_vec() {
         bdm_assert_debug ( Spoints.length() > idx - 1, "Out of support points" );
-        return Spoints ( ++idx );
+                idx++;
+                actvolume=volumes(idx);
+                actvec = Spoints ( idx );
+        return actvec;
     }
 
@@ -123 +162 @@
     void from_setting ( const Setting &set );
 
+        void validate(){
+                bdm_assert(volumes.length()==Spoints.length(),"dimension mismatch");
+                Npoints=Spoints.length();
+                if (Npoints>0)
+                        dim = Spoints(0).length();
+        }
     //! access function
     Array<vec> & _Spoints() {
@@ -134 +179 @@
 protected:
     //! grid - function support
-    rectangular_support sup;
+    shared_ptr<support_base> sup;
     //! function values on the grid
     vec values;
 public:
     //! constructor function
-    void set_support ( rectangular_support &sup0 ) {
+    void set_support ( shared_ptr<support_base> &sup0 ) {
         sup = sup0;
-        values = zeros ( sup.points() );
+        values = zeros ( sup->points() );
     }
     //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
@@ -151 +196 @@
     //! get value nearest to the given point
     double nearest_val ( const vec &val ) {
-        return values ( sup.linear_index ( sup.nearest_point ( val ) ) );
-    }
+                NOT_IMPLEMENTED(0.0);
+//        return values ( sup.linear_index ( sup.nearest_point ( val ) ) );
+    }
+
+        double integrate() ;
 
     vec eval ( const vec &val ) {