diff --git a/include/gridpp.h b/include/gridpp.h
index 4793e643..f60c249d 100644
--- a/include/gridpp.h
+++ b/include/gridpp.h
@@ -293,17 +293,34 @@ namespace gridpp {
             int max_points,
             bool allow_extrapolation=true);
 
-    vec2 optimal_interpolation_ensi(const Points& bpoints,
-            const vec2& background,       // T
-            // const vec2& other_background, // Precip
-            const vec2& original_background, // T
-            // const vec2& other_original_background,  // Precip
+    /** Optimal interpolation for an ensemble gridded field (alternative version)
+      * Work in progress
+      * @param bgrid Grid of background field
+      * @param background 3D vector of (left) background values to update (Y, X, E)
+      * @param background 3D vector of (LEFT) background values (Y, X, E) used to compute correlations
+      * @param obs_points observation points
+      * @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 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
+      * @param max_points Maximum number of observations to use inside localization zone; Use 0 to disable
+      * @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,
+            const vec3& background_l,
+            const vec3& background_L,
             const Points& obs_points,
-            const vec& obs,  // Precip
-            const vec& obs_standard_deviations, // Precip
-            const vec2& background_at_points, // Precip
-            const vec2& original_background_at_points, // Precip origin
+            const vec& obs,
+            const vec2& pbackground_r,
+            const vec2& pbackground_R,
             const StructureFunction& structure,
+            float var_ratios_or,
+            float std_ratios_lr,
+            float weigth,
             int max_points,
             bool allow_extrapolation=true);
 
diff --git a/src/api/oi_ensi.cpp b/src/api/oi_ensi.cpp
index 66f2551d..4f850088 100644
--- a/src/api/oi_ensi.cpp
+++ b/src/api/oi_ensi.cpp
@@ -112,8 +112,6 @@ vec3 gridpp::optimal_interpolation_ensi(const gridpp::Grid& bgrid,
 }
 vec2 gridpp::optimal_interpolation_ensi(const gridpp::Points& bpoints,
         const vec2& background,
-        const vec2& other_background,
-        const vec2& original_background,
         const gridpp::Points& points,
         const vec& pobs,  // gObs
         const vec& psigmas,   // pci
@@ -128,7 +126,6 @@ vec2 gridpp::optimal_interpolation_ensi(const gridpp::Points& bpoints,
     }
     if(background.size() != bpoints.size())
         throw std::invalid_argument("Input field is not the same size as the grid");
-    // TODO: Check new backgrounds
     if(pobs.size() != points.size())
         throw std::invalid_argument("Observations and points exception mismatch");
     if(psigmas.size() != points.size())
@@ -187,7 +184,6 @@ vec2 gridpp::optimal_interpolation_ensi(const gridpp::Points& bpoints,
     }
 
     // Calculate number of valid members
-    // TODO: Check the two other backgrounds
     int nValidEns = 0;
     ivec validEns;
     for(int e = 0; e < nEns; e++) {
@@ -452,7 +448,7 @@ vec2 gridpp::optimal_interpolation_ensi(const gridpp::Points& bpoints,
         int count = 0;
         for(int e = 0; e < nValidEns; e++) {
             int ei = validEns[e];
-            float value = original_background[y][ei];
+            float value = background[y][ei];
             if(gridpp::is_valid(value)) {
                 X(e) = value;
                 total += value;
diff --git a/src/api/oi_ensi_lr.cpp b/src/api/oi_ensi_lr.cpp
new file mode 100644
index 00000000..db5e5020
--- /dev/null
+++ b/src/api/oi_ensi_lr.cpp
@@ -0,0 +1,405 @@
+#include "gridpp.h"
+#include <math.h>
+#include <algorithm>
+#include <armadillo>
+#include <assert.h>
+#include <exception>
+
+using namespace gridpp;
+
+namespace {
+    typedef arma::mat mattype;
+    typedef arma::vec vectype;
+    typedef arma::cx_mat cxtype;
+
+    void check_vec(vec2 input, int Y, int X);
+    void check_vec(vec input, int S);
+
+    template<class T1, class T2> struct sort_pair_first {
+        bool operator()(const std::pair<T1,T2>&left, const std::pair<T1,T2>&right) {
+            return left.first < right.first;
+        };
+    };
+}
+
+template<class Matrix>
+void print_matrix(Matrix matrix) {
+       matrix.print(std::cout);
+}
+
+template void print_matrix< ::mattype>(::mattype matrix);
+template void print_matrix< ::cxtype>(::cxtype matrix);
+
+vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid,
+        const vec3& background_l, // left background: field to update
+        const vec3& background_L, // LEFT background: field used for dynamical correlations
+        const gridpp::Points& points,
+        const vec2& pobs,
+        const vec2& pbackground_r, // right background: new information at observations points
+        const vec2& pbackground_R, // RIGHT background: information at obs points for dynamical correlations
+        const gridpp::StructureFunction& structure,
+        float var_ratios_or, // eps2
+        float std_ratios_lr, // gamma
+        float weigth, // beta
+        int max_points,
+        bool allow_extrapolation) {
+    double s_time = gridpp::clock();
+
+    // Check input data
+    if(max_points < 0)
+        throw std::invalid_argument("max_points must be >= 0");
+
+    int nS = points.size();
+    if(nS == 0)
+        return background_l;
+
+    int nY = bgrid.size()[0];
+    int nX = bgrid.size()[1];
+
+    if(nY == 0 || nX == 0) {
+        std::stringstream ss;
+        ss << "Grid size (" << nY << "," << nX << ") cannot be zero";
+        throw std::invalid_argument(ss.str());
+    }
+
+    if(bgrid.get_coordinate_type() != points.get_coordinate_type()) {
+        throw std::invalid_argument("Both background grid and observations points must be of same coordinate type (lat/lon or x/y)");
+    }
+    // Check ensembles have consistent sizes
+    if(background_l.size() != nY || background_l[0].size() != nX) {
+        std::stringstream ss;
+        ss << "Input left field (" << background_l.size() << "," << background_l[0].size() << "," background_l[0][0].size() << ") is not the same size as the grid (" << nY << "," << nX << "," << nE << ")";
+        throw std::invalid_argument(ss.str());
+    }
+    int nE = background_l[0][0].size();
+    if(background_L.size() != nY || background_L[0].size() != nX || background_L[0][0].size() != nE) {
+        std::stringstream ss;
+        ss << "Input LEFT field (" << background_L.size() << "," << background_L[0].size() << "," background_L[0][0].size() << ") is not the same size as the grid (" << nY << "," << nX << "," << nE << ")";
+        throw std::invalid_argument(ss.str());
+    }
+    if(pbackground_r.size() != nS || pbackground_r[0].size() != nE) {
+        std::stringstream ss;
+        ss << "Input right field at observation location (" << pbackground_r.size() << "," << pbackground_r[0].size() << ") and points (" << nS "," << nE << ") size mismatch";
+        throw std::invalid_argument(ss.str());
+    }
+
+    if(pbackground_R.size() != nS || pbackground_R[0].size() != nE) {
+        std::stringstream ss;
+        ss << "Input RIGHT field at observation location (" << pbackground_R.size() << "," << pbackground_R[0].size() << ") and points (" << nS "," << nE << ") size mismatch";
+        throw std::invalid_argument(ss.str());
+    }
+    // Check observations have consistent size
+    if(pobs.size() != nS || pobs[0].size() != nE) {
+        std::stringstream ss;
+        ss << "Observations (" << pobs.size() << "," << pobs[0].size() << ") and points (" << nS "," << nE << ") size mismatch";
+        throw std::invalid_argument(ss.str());
+    }
+    
+    gridpp::Points bpoints = bgrid.to_points();
+    vec2 background_l1 = gridpp::init_vec2(nY * nX, nE);
+    vec2 background_L1 = gridpp::init_vec2(nY * nX, nE);
+    int count = 0;
+    for(int y = 0; y < nY; y++) {
+        for(int x = 0; x < nX; x++) {
+            for(int e = 0; e < nE; e++) {
+                background_l1[count][e] = background_l[y][x][e];
+                background_L1[count][e] = background_L[y][x][e];
+            }
+            count++;
+        }
+    }
+    vec2 output1 = optimal_interpolation_ensi_lr(bpoints, background_l1, background_L1, points, pobs, pbackground_r, pbackground_R, structure, var_ratios_or, std_ratios_lr, weigth, max_points, allow_extrapolation);
+    vec3 output = gridpp::init_vec3(nY, nX, nE);
+    count = 0;
+    for(int y = 0; y < nY; y++) {
+        for(int x = 0; x < nX; x++) {
+            for(int e = 0; e < nE; e++) {
+                output[y][x][e] = output1[count][e];
+            }
+            count++;
+        }
+    }
+    return output;
+}
+vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
+        const vec2& background_l, // left background: field to update
+        const vec2& background_L, // LEFT background: field used for dynamical correlations
+        const gridpp::Points& points,
+        const vec& pobs,  // gObs
+        const vec2& pbackground_r, // right background: new information at observations points
+        const vec2& pbackground_R, // RIGHT background: information at obs points for dynamical correlations
+        const gridpp::StructureFunction& structure,
+        float var_ratios_or, // eps2
+        float std_ratios_lr, // gamma
+        float weigth, // beta
+        int max_points,
+        bool allow_extrapolation) {
+    if(max_points < 0)
+        throw std::invalid_argument("max_points must be >= 0");
+    if(bpoints.get_coordinate_type() != points.get_coordinate_type()) {
+        throw std::invalid_argument("Both background and observations points must be of same coorindate type (lat/lon or x/y)");
+    }
+    if(background_l.size() != bpoints.size())
+        throw std::invalid_argument("Input left field is not the same size as the grid");
+    if(background_L.size() != bpoints.size())
+        throw std::invalid_argument("Input LEFT field is not the same size as the grid");
+    if(pobs.size() != points.size())
+        throw std::invalid_argument("Observations and points exception mismatch");
+    if(pbackground_r.size() != points.size())
+        throw std::invalid_argument("Background rigth and points size mismatch");
+    if(pbackground_R.size() != points.size())
+        throw std::invalid_argument("Background RIGTH and points size mismatch");
+
+    int nS = points.size();
+    if(nS == 0)
+        return background_l;
+
+    int mY = -1;  // Write debug information for this station index
+    bool diagnose = false;
+
+    int nY = background_l.size();
+    int nEns = background_l[0].size();
+
+    // Prepare output matrix
+    float missing_value = -99999.999;
+    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();
+
+    vec plats = points.get_lats();
+    vec plons = points.get_lons();
+    vec pelevs = points.get_elevs();
+    vec plafs = points.get_lafs();
+
+
+    // Create all objects of type Point (to save time on creation later)
+    std::vector<Point> point_vec;
+    point_vec.reserve(nS);
+    for(int i = 0; i < nS; i++) {
+        point_vec.push_back(points.get_point(i));
+    }
+
+    // Calculate number of valid members
+    // NOTE: all backgrounds must be simultaneously valid
+    int nValidEns = 0;
+    ivec validEns;
+    for(int e = 0; e < nEns; e++) {
+        int numInvalid = 0;
+        for(int y = 0; y < nY; y++) {
+            float value_l = background_l[y][e];
+            float value_L = background_L[y][e];
+            if(!gridpp::is_valid(value_l) || !gridpp::is_valid(value_L))
+                numInvalid++;
+        }
+        for(int i = 0; i < nS; i++) {
+            float value_r = pbackground_r[y][e];
+            float value_R = pbackground_R[y][e];
+            if(!gridpp::is_valid(value_r) || !gridpp::is_valid(value_R))
+                numInvalid++;
+        }
+        if(numInvalid == 0) {
+            validEns.push_back(e);
+            nValidEns++;
+        }
+    }
+    if(nValidEns == 0)
+        return background_l;
+
+    // gZ_R(nY, nValidEns): used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observationsCompute Y
+    vec2 gZ_R = gridpp::init_vec2(nY, nValidEns); // useful to compute dynamical correlations
+    for(int i = 0; i < nS; i++) {
+        vec pbackgroundValid_R(nValidEns);
+        for(int e = 0; e < nValidEns; e++) {
+            int ei = validEns[e];
+            pbackgroundValid_R[e] = pbackground_R[i][ei];
+        }
+        float mean = gridpp::calc_statistic(pbackgroundValid_R, gridpp::Mean);
+        float std = gridpp::calc_statistic(pbackgroundValid_R, gridpp::Std);
+        if(!gridpp::is_valid(mean) || !gridpp::is_valid(std)) {
+            for(int e = 0; e < nValidEns; e++) 
+                gZ_R[i][e] = 0;
+        }
+        else {
+            if(std == 0) {
+                for(int e = 0; e < nValidEns; e++) 
+                    gZ_R[i][e] = 0;
+            }
+            else {
+                for(int e = 0; e < nValidEns; e++) 
+                    gZ_R[i][e] = 1 / sqrt(nValidEns-1) * (pbackgroundValid_R[e] - mean) / std;
+            }
+        }
+    } // end loop over all observations
+
+    // This causes segmentation fault when building the gridpp pypi package
+    // 1) Tested removing num_condition_warning and num_real_part_warning from the parallel loop
+    //    but it doesnt seem to help
+    // #pragma omp parallel for
+    for(int y = 0; y < nY; y++) {
+        float lat = blats[y];
+        float lon = blons[y];
+        float elev = belevs[y];
+        float laf = blafs[y];
+        Point p1 = bpoints.get_point(y);
+        float localizationRadius = structure.localization_distance(p1);
+
+        // Create list of locations for this gridpoint
+        ivec lLocIndices0 = points.get_neighbours(lat, lon, localizationRadius);
+        if(lLocIndices0.size() == 0) {
+            // If we have too few observations though, then use the background
+            continue;
+        }
+        ivec lLocIndices;
+        lLocIndices.reserve(lLocIndices0.size());
+        std::vector<std::pair<float,int> > lRhos0;
+
+        // Calculate gridpoint to observation rhos
+        lRhos0.reserve(lLocIndices0.size());
+        std::vector<Point> p2;
+        p2.reserve(lLocIndices0.size());
+        for(int i = 0; i < lLocIndices0.size(); i++) {
+            int index = lLocIndices0[i];
+            p2.push_back(point_vec[index]);
+        }
+        vec rhos = structure.corr_background(p1, p2);
+        for(int i = 0; i < lLocIndices0.size(); i++) {
+            int index = lLocIndices0[i];
+            if(gridpp::is_valid(pobs[index])) {
+                if(rhos_static_lr[i] > 0) {
+                    lRhos0.push_back(std::pair<float,int>(rhos_static_lr[i], i));
+                }
+            }
+        }
+
+        // Make sure we don't use too many observations
+        arma::vec lRhos;
+        if(max_points > 0 && lRhos0.size() > max_points) {
+            // If sorting is enabled and we have too many locations, then only keep the best ones based on rho.
+            // Otherwise, just use the last locations added
+            lRhos = arma::vec(max_points);
+            std::sort(lRhos0.begin(), lRhos0.end(), ::sort_pair_first<float,int>());
+            for(int i = 0; i < max_points; i++) {
+                // The best values start at the end of the array
+                int index = lRhos0[lRhos0.size() - 1 - i].second;
+                lLocIndices.push_back(lLocIndices0[index]);
+                lRhos(i) = lRhos0[lRhos0.size() - 1 - i].first;
+            }
+        }
+        else {
+            lRhos = arma::vec(lRhos0.size());
+            for(int i = 0; i < lRhos0.size(); i++) {
+                int index = lRhos0[i].second;
+                lLocIndices.push_back(lLocIndices0[index]);
+                lRhos(i) = lRhos0[i].first;
+            }
+        }
+
+        int lS = lLocIndices.size();
+        if(lS == 0) {
+            // If we have too few observations though, then use the background
+            continue;
+        }
+
+        vectype lElevs(lS);
+        vectype lLafs(lS);
+        for(int i = 0; i < lLocIndices.size(); i++) {
+            int index = lLocIndices[i];
+            lElevs[i] = pelevs[index];
+            lLafs[i] = plafs[index];
+        }
+
+        // lX_L: used to compute ensemble-based background correlations between yth gridpoint and observations
+        mattype lX_L(1, nValidEns, arma::fill::zeros);
+        vec backgroundValid_L(nValidEns);
+        for(int e = 0; e < nValidEns; e++) {
+            int ei = validEns[e];
+            backgroundValid_L[e] = background_L[y][ei];
+        }
+        float mean = gridpp::calc_statistic(backgroundValid_L, gridpp::Mean);
+        float std = gridpp::calc_statistic(backgroundValid_L, gridpp::Std);
+        if(gridpp::is_valid(mean) && gridpp::is_valid(std) && std != 0) {
+            for(int e = 0; e < nValidEns; e++) 
+                lX_L[0][e] = 1 / sqrt(nValidEns-1) * (backgroundValid_L[e] - mean) / std;
+        }
+        // lZ_R: used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observations
+        mattype lZ_R(lS, nValidEns);
+        // Innovation: Observation - Background_r
+        mattype lInnov(lS, nValidEns, arma::fill::zeros);
+        // lR_dd: Observation error correlation matrix
+        mattype lR_dd(lS, lS, arma::fill::zeros);
+        // lLoc1D: localization for ensemble-based background correlations between yth gridpoint and observations
+        mattype lLoc1D(1, lS, arma::fill::zeros); 
+        // lLoc2D: localization for ensemble-based background correlations among observations
+        mattype lLoc2D(lS, lS, arma::fill::zeros); 
+        for(int i = 0; i < lS; i++) {
+            // lR_dd: diagonal observation error correlation matrix
+            lR_dd(i, i) = var_ratios_or;
+            // lLoc1D between yth gridpoint and ith observation computed before
+            lLoc1D(0, i) = lRhos(i);
+            // compute lZ_R and lInnov
+            int index = lLocIndices[i];
+            for(int e = 0; e < nValidEns; e++) {
+                lZ_R(i, e) = gZ_R[index][e];
+                int ei = validEns[e];
+                lInnov(i, ei) = pobs[index][ei] - pbackground_r[index][ei];
+            }
+            // compute lLoc2D
+            Point p1 = point_vec[index];
+            std::vector<Point> p2(lS, Point(0, 0));
+            for(int j = 0; j < lS; j++) {
+                int index_j = lLocIndices[j];
+                p2[j] = point_vec[index_j];
+            }
+            vec corr = structure.corr(p1, p2);
+            for(int j = 0; j < lS; j++) 
+                lLoc2D(i, j) = corr[j];
+        } // end loop over closer observations
+        // lr_lr(1, lS): ensemble-based background correlations between yth gridpoint and observations
+        //  note: correlation between the LEFT background and the RIGHT background (..._lr)
+        mattype lr_lr = lLoc1D % (lX_L * lZ_R.t());
+        // lR_rr(lS, lS): ensemble-based background correlations among observations
+        //  note: correlation between the RIGTH background and the RIGHT background (..._rr)
+        mattype lR_rr = lLoc2D % (lZ_R * lZ_R.t());
+        // 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);
+
+        ///////////////////////////////
+        // Anti-extrapolation filter //
+        ///////////////////////////////
+        if(!allow_extrapolation) {
+            vectype MaxIncEns = arma::max(lInnov);
+            vectype MinIncEns = arma::min(lInnov);
+
+            for(int e = 0; e < nValidEns; e++) {
+                float increment = dx[e];
+                float MaxInc = MaxIncEns[e];
+                float MinInc = MinIncEns[e];
+                if(maxInc > 0 && increment > maxInc) {
+                   increment = maxInc;
+                }
+                else if(maxInc < 0 && increment > 0) {
+                   increment = maxInc;
+                }
+                else if(minInc < 0 && increment < minInc) {
+                   increment = minInc;
+                }
+                else if(minInc > 0 && increment < 0) {
+                   increment = minInc;
+                }
+                dx[e] = increment;
+            }
+        }
+
+        // Compute the analysis as the updated background_l
+        for(int e = 0; e < nValidEns; e++) {
+            int ei = validEns[e];
+            output[y][ei] = background_l[y][ei] + dx[e];
+        }
+    } // end loop over gridpoint 
+    return output;
+}