modified such that it works with the R bindings

Author: Cristian Lussana

include/gridpp.h

@@ -310,7 +310,7 @@ namespace gridpp {
       * @param allow_extrapolation Allow OI to extrapolate increments outside increments at observations
       * @returns 3D vector of analised values (Y, X, E)
     */
-    vec3 optimal_interpolation_ensi_lr(const Grid& bgrid,
+/*    vec3 optimal_interpolation_ensi_lr(const Grid& bgrid,
             const vec3& background_l,
             const vec3& background_L,
             const Points& obs_points,
@@ -322,7 +322,7 @@ namespace gridpp {
             float std_ratios_lr,
             float weigth,
             int max_points,
-            bool allow_extrapolation=true);
+            bool allow_extrapolation=true); */
 
     /** Optimal interpolation for an ensemble gridded field (alternative version)
       * Work in progress
@@ -333,7 +333,9 @@ namespace gridpp {
       * @param obs 2D vector of perturbed observations (S, E)
       * @param background 3D vector of (right) background values used to compute innovations (Y, X, E)
       * @param background 3D vector of (RIGHT) background values (Y, X, E) used to compute correlations
-      * @param structure Structure function
+      * @param dh structure Structure function par 1
+      * @param dz structure Structure function par 2
+      * @param dw structure Structure function par 3
       * @param variance_ratio (ratio of observation to right background error variance)
       * @param standard deviation ratio (ratio of left to right background error standard deviation)
       * @param weight given to the analysis increment
@@ -348,7 +350,10 @@ namespace gridpp {
             const vec2& obs,
             const vec2& pbackground_r,
             const vec2& pbackground_R,
-            const StructureFunction& structure,
+/*            const StructureFunction& structure, */
+            float dh,
+            float dz,
+            float dw,
             float var_ratios_or,
             float std_ratios_lr,
             float weigth,
@@ -1591,6 +1596,8 @@ namespace gridpp {
     float test_vec2_argout(vec2& distances);
 
     void test_not_implemented_exception();
+/*    vec2 test_args_for_R(const Points& bpoints, const StructureFunction& structure, const vec2& background);*/
+    vec2 test_args_for_R(const Points& bpoints, const vec2& background);
 
     /** Default value used to fill array in SWIG testing functions. Not useful for any other purpose. */
     static const float swig_default_value = -1;

src/api/oi_ensi_lr.cpp

@@ -30,7 +30,7 @@ void print_matrix(Matrix matrix) {
 template void print_matrix< ::mattype>(::mattype matrix);
 template void print_matrix< ::cxtype>(::cxtype matrix);
 
-vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid,
+/*vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid,
         const vec3& background_l, 
         const vec3& background_L, 
         const gridpp::Points& points,
@@ -120,15 +120,18 @@ vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid,
         }
     }
     return output;
-}
+} */
 vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
         const vec2& background_l, 
         const vec2& background_L, 
         const gridpp::Points& points,
         const vec2& pobs,  
         const vec2& pbackground_r, 
         const vec2& pbackground_R,
-        const gridpp::StructureFunction& structure,
+/*        const gridpp::StructureFunction& structure, it creates problem with R bindings */
+        float dh,
+        float dz,
+        float dw,
         float var_ratios_or,
         float std_ratios_lr,
         float weigth,
@@ -150,6 +153,9 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
     if(pbackground_R.size() != points.size())
         throw std::invalid_argument("Background RIGTH and points size mismatch");
 
+    float hmax = 7 * dh;
+    BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
+
     int nS = points.size();
     if(nS == 0)
         return background_l;
@@ -162,7 +168,8 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
 
     // Prepare output matrix
     float missing_value = -99999.999;
-    vec2 output = gridpp::init_vec2(nY, nEns, missing_value);
+/*    vec2 output = gridpp::init_vec2(nY, nEns, missing_value); */
+    vec2 output = background_l;
 
     vec blats = bpoints.get_lats();
     vec blons = bpoints.get_lons();
@@ -366,8 +373,7 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
         // lK(1, lS): Kalman gain
         mattype lK = lr_lr * arma::inv(lR_rr + lR_dd);
         // dx(1, nValidEns): analysis increment 
-        vectype dx = std_ratios_lr * weigth * (lK * lInnov);
-
+        mattype dx = std_ratios_lr * weigth * (lK * lInnov);
         ///////////////////////////////
         // Anti-extrapolation filter //
         ///////////////////////////////
@@ -400,6 +406,17 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
             int ei = validEns[e];
             output[y][ei] = background_l[y][ei] + dx[e];
         }
+        // debug
+/*        for(int i = 0; i < lS; i++) {
+            // compute lZ_R and lInnov
+            int index = lLocIndices[i];
+            std::cout << i << " backg_r obs " << pbackground_r[index][0] << " " << pobs[index][0] << std::endl;
+        } // end loop over closer observations
+        for(int e = 0; e < nValidEns; e++) {
+            int ei = validEns[e];
+            std::cout << ei << " backg_l analysis " << background_l[y][ei] << " " << output[y][ei] << std::endl;
+        } */
+
     } // end loop over gridpoint 
     return output;
 }

src/api/swig.cpp

@@ -98,3 +98,36 @@ float gridpp::test_vec2_argout(vec2& distances) {
 void gridpp::test_not_implemented_exception() {
     throw gridpp::not_implemented_exception();
 }
+vec2 gridpp::test_args_for_R(const gridpp::Points& bpoints,
+/*                             const gridpp::StructureFunction& structure,*/
+                             const vec2& background) {
+    int nS = bpoints.size();
+    int nEns = background[0].size();
+    int nY = background.size();
+    float missing_value = -99999.999;
+/*    gridpp::StructureFunction structure = gridpp::BarnesStructure(10000, 0, 0, 0);*/
+    BarnesStructure structure = BarnesStructure(10000,1000000,0.5,100000);
+    vec2 output = gridpp::init_vec2(nY, nEns, missing_value);
+    vec blats = bpoints.get_lats();
+    vec blons = bpoints.get_lons();
+    vec belevs = bpoints.get_elevs();
+    vec blafs = bpoints.get_lafs();
+    Point p2 = bpoints.get_point(0);
+    for(int i = 0; i < nS; i++) {
+        Point p1 = bpoints.get_point(i);
+        float localizationRadius = structure.localization_distance(p1);
+        float rhos = structure.corr_background(p1, p2);
+        std::cout << "----------------------------" << std::endl;
+        std::cout << i << " localizationRadius " << localizationRadius << std::endl;
+        std::cout << i << " rhos " << rhos << std::endl;
+        std::cout << i << " blats " << blats[i] << std::endl;
+        std::cout << i << " blons " << blons[i] << std::endl;
+        std::cout << i << " belevs " << belevs[i] << std::endl;
+        std::cout << i << " blafs " << blafs[i] << std::endl;
+        for(int j = 0; j < background[i].size(); j++) {
+            std::cout << j << " bkg " << background[i][j] << std::endl;
+            output[i][j] = background[i][j];
+        }
+    }
+    return output;
+}