Repository: metno/gridpp (language: python)

Issue:
Implement calc_gradient
Implement min/max gradient function for now. Later include linear regression.

Relevant files: input/gridpp.h and src/api/calc_gradient.cpp
- [x]  Remove grid as an input to function
- [x] Rename radius to halfwidth
- [x]  Implement logic in calc_gradient.cpp (look at neighbourhood.cpp for an example)
- [x]  Implement tests in tests/calc_gradient_test.py (see tests/wind_test.py for an example)


Discussion / hints:

Fix patch:
diff --git a/extras/SWIG/gridpp.i b/extras/SWIG/gridpp.i
index 49bdd382..210692e6 100644
--- a/extras/SWIG/gridpp.i
+++ b/extras/SWIG/gridpp.i
@@ -53,6 +53,7 @@
 %apply int& OUTPUT { int& Y1_out };
 %apply int& OUTPUT { int& X2_out };
 %apply int& OUTPUT { int& Y2_out };
+%apply std::vector<std::vector<std::vector<float> > >& OUTPUT { std::vector<std::vector<std::vector<float> > >& output };
 
 %{
 #include "gridpp.h"
diff --git a/include/gridpp.h b/include/gridpp.h
index 0405867e..60c9c33a 100644
--- a/include/gridpp.h
+++ b/include/gridpp.h
@@ -97,6 +97,13 @@ namespace gridpp {
         Cartesian = 1,     /**< X and Y */
     };
 
+    /** Types of methods to calculate the gradient*/
+    enum GradientType {
+        MinMax = 0,
+        LinearRegression = 10,
+    };
+
+
     /** **************************************
      * @name Data assimilation methods
      * Functions that merge observations with a background field
@@ -768,7 +775,20 @@ namespace gridpp {
      *  @param min_range Minimum range of base to compute gradient
      *  @param default_gradient Use this gradient if minimum number is not met
     */
-    vec2 calc_gradient(const Grid& grid, const vec2& base, const vec2& values, int radius, int min_num=2, float min_range=0, float default_gradient=0);
+    vec2 calc_gradient(const vec2& base, const vec2& values, int halfwidth, int min_num=2, float min_range=0, float default_gradient=0);
+
+    /** Compute downscale 
+    *@param igrid input grid
+    *@param ogrid output grid 
+    *@param ivalues temperature values from igrid
+    *@param elev_gradient elevation gradient
+    *@param laf_gradient land area fraction gradient
+    */
+    vec2 gradient(const Grid& igrid, const Grid& ogrid, const vec2& ivalues,  const vec2& elev_gradient, const vec2& laf_gradient);
+    
+    /* Debug method. Delete later
+    */
+    vec3 correction(const Grid& igrid, const Grid& ogrid, const vec2& ivalues,  const vec2& elev_gradient, const vec2& laf_gradient);
 
     /** Check if the grid is the same size as the 2D vector */
     bool compatible_size(const Grid& grid, const vec2& v);
diff --git a/src/api/.calc_gradient.cpp.swo b/src/api/.calc_gradient.cpp.swo
new file mode 100644
index 00000000..0bed262c
Binary files /dev/null and b/src/api/.calc_gradient.cpp.swo differ
diff --git a/src/api/calc_gradient.cpp b/src/api/calc_gradient.cpp
index 1e03e639..19bca13c 100644
--- a/src/api/calc_gradient.cpp
+++ b/src/api/calc_gradient.cpp
@@ -1,7 +1,150 @@
 #include "gridpp.h"
+#include <iostream>
 
 using namespace gridpp;
 
-vec2 gridpp::calc_gradient(const gridpp::Grid& grid, const vec2& base, const vec2& values, int radius, int min_num, float min_range, float default_gradient) {
-    throw gridpp::not_implemented_exception();
-}
+/*GradientType{
+*   minMax = 0;
+*   LinearRegression = 10;
+};*/
+
+vec2 gridpp::calc_gradient(const vec2& base, const vec2& values, //GradientType gradientType,
+    int halfwidth , int min_num, float min_range, float default_gradient){
+
+    if(halfwidth < 0)
+        throw std::invalid_argument("halfwidth must be positive");
+
+    //std::cout << " Starting ";
+
+    int gradientType = 0;
+
+    //if(halfwidth == 0)
+    //    throw std::invalid_argument("Halwidth cannot be 0; must be positive integer");
+    //if(base.size() == 0)
+    //    throw std::invalid_argument("Base input has no size");
+
+
+    int nY = base.size();
+    int nX = base[0].size();
+
+    vec2 output = gridpp::init_vec2(base.size(), base[0].size());
+
+    // Min Max Gradient 
+    if(gradientType == 0){
+        for(int y = 0; y  < base.size(); y++) {
+            for(int x = 0; x < base[y].size(); x++){
+                bool start = true;
+
+                float current_max = 0;
+                float current_min = 0;
+                int I_maxBase_Y = 0;
+                int I_maxBase_X = 0;
+                int I_minBase_Y = 0;
+                int I_minBase_X = 0;
+                for(int yy = std::max(0, y - halfwidth); yy <= std::min(nY - 1, y + halfwidth); yy++){               
+                    for(int xx = std::max(0, x - halfwidth); xx <= std::min(nX - 1, x + halfwidth); xx++){
+                        float current_base = base[yy][xx];
+                        if(!gridpp::is_valid(current_base)){
+                            continue; 
+                        }
+
+                        if(start){
+                            start = false;
+                            current_max = current_base;
+                            I_maxBase_Y = yy;
+                            I_maxBase_X = xx;
+                            current_min = current_base;
+                            I_minBase_Y = yy;
+                            I_minBase_X = xx;
+                        }
+
+                        else if(current_base > current_max){
+                            current_max = current_base;
+                            I_maxBase_Y = yy;
+                            I_maxBase_X = xx;
+                        }
+                        else if(current_base < current_min){
+                            current_min = current_base;
+                            I_minBase_Y = yy;
+                            I_minBase_X = xx;
+                        }
+                    }       
+                }
+                if(!gridpp::is_valid(current_max) || !gridpp::is_valid(current_min)){
+                    output[y][x] = default_gradient;
+                }
+
+                else if(abs(current_max - current_min) < min_range){
+                    output[y][x] = default_gradient;
+                }
+                else{
+                    float diffBase = current_max - current_min;
+                    float diffValues = values[I_maxBase_Y][I_maxBase_X] - values[I_minBase_Y][I_minBase_X];
+                    output[y][x] = diffValues / diffBase;
+                }
+            }
+        }
+    }
+    // Linear Regression Gradient
+    else if(gradientType == 10){
+        for(int y = 0; y  < base.size(); y++) {
+            for(int x = 0; x < base[y].size(); x++){             
+                bool start = true;
+                float current_max = 0;
+                float current_min = 0;
+
+                float meanX = 0; 
+                float meanY = 0;
+                float meanXY = 0;
+                float meanXX = 0;
+
+                int counter = 0;
+                
+                for(int yy = std::max(0, y - halfwidth); yy <= std::min(nY - 1, y + halfwidth); yy++){               
+                    for(int xx = std::max(0, x - halfwidth); xx <= std::min(nX - 1, x + halfwidth); xx++){
+                        float current_base = base[yy][xx];
+                        float current_values = values[yy][xx];
+
+                        if(start){
+                            start = false;
+                            current_max = current_base;
+                            current_min = current_base;
+                        }
+
+                        else if(current_base > current_max){
+                            current_max = current_base;
+                        }
+
+                        else if(current_base < current_min){
+                            current_min = current_base;
+                        }
+
+                        meanX = meanX + current_base;
+                        meanY = meanY + current_values;
+                        meanXX = meanXX + current_base*current_base;
+                        meanXY = meanXY + current_base*current_values;
+                        counter = counter + 1;
+                    }       
+                }
+
+                if(abs(current_max - current_min) < min_range){
+                    output[y][x] = default_gradient;
+                }
+
+                else{
+                    meanX = meanX / counter;
+                    meanY = meanY / counter;
+                    meanXX = meanXX / counter;
+                    meanXY = meanXY / counter;
+                    if(meanXX - meanX*meanX != 0){
+                        output[y][x] = (meanXY - meanX*meanY)/(meanXX - meanX*meanX);
+                    }
+                    else{
+                        output[y][x] = default_gradient;
+                    }
+                }
+            }
+        }        
+    }
+    return output;
+}
\ No newline at end of file
diff --git a/src/api/gradient.cpp b/src/api/gradient.cpp
new file mode 100644
index 00000000..ccf52f91
--- /dev/null
+++ b/src/api/gradient.cpp
@@ -0,0 +1,128 @@
+#include "gridpp.h"
+
+using namespace gridpp;
+
+
+
+vec2 gridpp::gradient(const Grid& igrid, const Grid& ogrid, const vec2& ivalues,  const vec2& elev_gradient, const vec2& laf_gradient){
+    //Inputs
+    //vec2 ivalues; 
+    vec2 ilats = igrid.get_lats();
+    vec2 ilons = igrid.get_lons();
+    vec2 ielevs = igrid.get_elevs();
+    vec2 ilafs = igrid.get_lafs();
+
+    vec2 i_lafs_gradient = laf_gradient;
+    vec2 i_elevs_gradient = elev_gradient;
+
+    //Outputs
+    //vec2 o_values = ogrid::nearest(igrid, ogrid, ivalues); Dont need this?
+    vec2 olats = ogrid.get_lats();
+    vec2 olons = ogrid.get_lons();
+    vec2 oelevs = ogrid.get_elevs();
+    vec2 olafs = ogrid.get_lafs();
+
+
+
+    // Sizes;
+    //int nY = ogrid.size();
+    //int nX = ogrid[0].size(); 
+
+    int nY = ogrid.size()[0]; 
+    int nX = ogrid.size()[1];
+
+    vec2 output = gridpp::init_vec2(nY, nX);
+
+    for(int y = 0; y < nY ; y++){
+        for(int x = 0; x < nX; x++){
+            //collect index for nearest neighbour
+            ivec indices = igrid.get_nearest_neighbour(olats[y][x], olons[y][x]);
+
+            //Collect output and input LAF and elevations
+            float olaf = olafs[y][x];
+            float ilaf = ilafs[indices[0]][indices[1]];
+            float oelev = oelevs[y][x];
+            float ielev = ielevs[indices[0]][indices[1]];
+
+            //Calculate LAF and elevation difference between output and input
+            float laf_diff = olaf - ilaf; 
+            float elev_diff = oelev - ielev;
+
+            float laf_correction = i_lafs_gradient[indices[0]][indices[1]]*laf_diff;
+            float elev_correction = i_elevs_gradient[indices[0]][indices[1]]*elev_diff;
+
+            //Calculate temperature
+            float temp = ivalues[indices[0]][indices[1]] + laf_correction + elev_correction;
+
+            //Assign value to output
+            output[y][x] = temp;
+        }
+    }
+    return output;
+}
+
+
+vec3 gridpp::correction(const Grid& igrid, const Grid& ogrid, const vec2& ivalues,  const vec2& elev_gradient, const vec2& laf_gradient){
+    //Inputs
+    //vec2 ivalues; 
+
+    vec2 ilats = igrid.get_lats();
+    vec2 ilons = igrid.get_lons();
+    vec2 ielevs = igrid.get_elevs();
+    vec2 ilafs = igrid.get_lafs();
+
+    vec2 i_lafs_gradient = laf_gradient;
+    vec2 i_elevs_gradient = elev_gradient;
+
+    //Outputs
+    //vec2 o_values = ogrid::nearest(igrid, ogrid, ivalues); Dont need this?
+    vec2 olats = ogrid.get_lats();
+    vec2 olons = ogrid.get_lons();
+    vec2 oelevs = ogrid.get_elevs();
+    vec2 olafs = ogrid.get_lafs();
+
+    // Sizes;
+    int nY = ogrid.size()[0];
+    int nX = ogrid.size()[1];
+
+    vec3 output = gridpp::init_vec3(nY, nX, 9);
+
+    for(int y = 0; y < nY ; y++){
+        for(int x = 0; x < nX; x++){
+            //collect index for nearest neighbour
+            ivec indices = igrid.get_nearest_neighbour(olats[y][x], olons[y][x]);
+
+            //Collect output and input LAF and elevations
+            float olaf = olafs[y][x];
+            float ilaf = ilafs[indices[0]][indices[1]];
+            //ilaf = ilaf/100.0;
+            float oelev = oelevs[y][x];
+            float ielev = ielevs[indices[0]][indices[1]];
+
+            //Calculate LAF and elevation difference between output and input
+            float laf_diff = olaf - ilaf;
+            float elev_diff = oelev - ielev;
+            
+            float laf_correction = i_lafs_gradient[indices[0]][indices[1]]*laf_diff;
+            float elev_correction = i_elevs_gradient[indices[0]][indices[1]]*elev_diff;
+
+            //Calculate temperature
+            float temp = ivalues[indices[0]][indices[1]] + laf_correction + elev_correction;
+            float temp_laf_c = ivalues[indices[0]][indices[1]] + laf_correction;
+            float temp_elev_c = ivalues[indices[0]][indices[1]]+ elev_correction;
+
+
+            //Assign value to output
+            output[y][x][0] = temp;
+            output[y][x][1] = laf_diff;
+            output[y][x][2] = elev_diff;
+            output[y][x][3] = laf_correction;
+            output[y][x][4] = elev_correction;
+            output[y][x][5] = olaf;
+            output[y][x][6] = ilaf;
+            output[y][x][7] = temp_laf_c;
+            output[y][x][8] = temp_elev_c;                        
+        }
+    }
+    return output;
+}
\ No newline at end of file


Test patch:
diff --git a/tests/benchmark.py b/tests/benchmark.py
index 305851dc..6701a894 100644
--- a/tests/benchmark.py
+++ b/tests/benchmark.py
@@ -7,7 +7,7 @@
 def main():
     parser = argparse.ArgumentParser(description='Runs gridpp benchmarks for processing performance')
     parser.add_argument('-j', type=int, help='Run multiple cores', dest='num_cores')
-    parser.add_argument('-s', type=int, default=1, help='Enlarge the inputs by this scaling factor to run a bigger test', dest='scaling')
+    parser.add_argument('-s', type=float, default=1, help='Enlarge the inputs by this scaling factor to run a bigger test', dest='scaling')
     parser.add_argument('-n', type=int, default=1, help='Number of iterations to average over', dest='iterations')
     parser.add_argument('-t', help='Run only this function', dest="function")
 
@@ -23,9 +23,9 @@ def main():
     np.random.seed(1000)
 
     for i in [10, 50, 100, 200, 500, 1000, 2000, 10000]:
-        input[i] = np.random.rand(i * args.scaling, i)*10
+        input[i] = np.random.rand(int(i * args.scaling), i)*10
     for i in [10, 50, 100, 200, 500, 1000]:
-        grids[i] = gridpp.Grid(*np.meshgrid(np.linspace(0, 1, i), np.linspace(0, 1, i * args.scaling)))
+        grids[i] = gridpp.Grid(*np.meshgrid(np.linspace(0, 1, i), np.linspace(0, 1, int(i * args.scaling))))
     for i in [1000]:
         # points[i] = gridpp.Points(np.linspace(0, 1, i), np.zeros(i))
         points[i] = gridpp.Points(np.random.rand(i) * 10, np.random.rand(i) * 10)
@@ -34,7 +34,7 @@ def main():
     quantile = 0.5
     thresholds = np.linspace(0, 1, 11)
     run = collections.OrderedDict()
-    run[(gridpp.Grid, "1000²")] = {"expected": 0.74, "args":np.meshgrid(np.linspace(0, 1, 1000 * args.scaling), np.linspace(0, 1, 1000))}
+    run[(gridpp.Grid, "1000²")] = {"expected": 0.74, "args":np.meshgrid(np.linspace(0, 1, int(1000 * args.scaling)), np.linspace(0, 1, 1000))}
     run[(gridpp.neighbourhood, "10000²")] = {"expected": 2.05, "args":(np.zeros([10000, 10000]), radius, gridpp.Mean)}
     run[(gridpp.neighbourhood,"2000² max")] = {"expected": 0.99, "args":(input[2000], radius, gridpp.Max)}
     run[(gridpp.neighbourhood_quantile_fast, "2000²")] = {"expected": 1.23, "args":(input[2000], quantile, radius, thresholds)}
@@ -46,6 +46,8 @@ def main():
     run[(gridpp.optimal_interpolation, "1000² 1000")] = {"expected": 1.57, "args":(grids[1000],
         input[1000], points[1000], np.zeros(1000), np.ones(1000), np.ones(1000), structure, 20)}
     run[(gridpp.dewpoint, "1e7")] = {"expected": 0.53, "args":(np.zeros(10000000) + 273.15, np.zeros(10000000))}
+    run[(gridpp.gradient, "1000²")] = {"expected": 3.55, "args": (grids[1000], grids[1000], np.zeros([1000,1000]), np.zeros([1000,1000]), np.zeros([1000,1000]))}
+    run[(gridpp.calc_gradient, "1000²")] = {"expected": 1.55, "args": (np.zeros([1000,1000]), np.zeros([1000,1000]), 3, 0, 0, 0)}
 
     print("Gridpp version %s" % gridpp.version())
     if args.num_cores is not None:
@@ -89,8 +91,8 @@ def main():
         for num_core in num_cores:
             timings[num_core] /= args.iterations
 
-        diff = (timings[1] - run[key]["expected"] * args.scaling) / (run[key]["expected"]  * args.scaling) * 100
-        string = "%-36s %8.2f %8.2f %8.2f %%" % (name, run[key]["expected"] * args.scaling, timings[1], diff)
+        diff = (timings[1] - int(run[key]["expected"] * args.scaling)) / int((run[key]["expected"]  * args.scaling) * 100)
+        string = "%-36s %8.2f %8.2f %8.2f %%" % (name, int(run[key]["expected"] * args.scaling), timings[1], diff)
         if args.num_cores is not None:
             scaling = timings[1] / timings[args.num_cores] / args.num_cores
             expected = timings[1] / args.num_cores
diff --git a/tests/gradient_test.py b/tests/gradient_test.py
new file mode 100644
index 00000000..a5c97e6d
--- /dev/null
+++ b/tests/gradient_test.py
@@ -0,0 +1,76 @@
+import gridpp
+import unittest
+import numpy as np 
+import random
+
+class Test(unittest.TestCase):	
+    def setUp(self):
+        self.input_base = [[0 ,0 , 0.2, 0.4], [0, 0.2, 0.4, 0.6], [0.2, 0.4,0.6, 0.8],[0.4, 0.6, 0.8, 1.0]]
+        #self.input_values = [[11, 11, 12, 14], [11, 12, 14, 18],[12, 14, 18, 26], [14, 18, 26, 42]] 
+        self.input_values = [[10,10,12,14], [10,12,14,16],[12,14,16,18],[14,16,18,20]] 
+
+    def test_size_of_inputs(self):
+        """ Checks that the input size of input and output gradient matrix are at the same size""" 
+        np.testing.assert_array_equal(len(self.input_base), len( gridpp.calc_gradient(self.input_base, self.input_values, 3))) 
+
+        #np.testing.assert_array_equal(self.input_base.shape, gridpp.calc_gradient(self.input_base, self.input_values, 3).shape)
+        #np.testing.assert_array_equal(self.input_values.shape, gridpp.calc_gradient(self.input_base, self.input_values, 3).shape)
+
+    def test_gradient(self):
+        """ tests that inputs give the correct gradient"""
+        self.input_values = [[10,10,12,14], [10,12,14,16],[12,14,16,18],[14,16,18,20]]
+        output = gridpp.calc_gradient(self.input_base, self.input_values, 3)
+        np.testing.assert_array_almost_equal(output, np.ones(output.shape)*10, decimal=3)
+
+
+    def test_nan_input(self):
+        """ Check gradient of inputs, where bases and values are nan, output should also be nan like base and values"""
+        self.input_base =  [[np.nan, np.nan, np.nan, np.nan]]*4
+        np.testing.assert_array_almost_equal(gridpp.calc_gradient(self.input_base,self.input_values, 3), np.nan)
+
+    def test_nan_ouput(self):
+        """ If an matrix contains a field with no changing base numberrs, test the output is nan"""
+        self.input_base = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
+        output = gridpp.calc_gradient(self.input_base, self.input_values, 3)
+        np.testing.assert_equal(output[0,0], np.nan)       
+
+    def test_empty_input(self):
+        """ Check empty input returns and empty output. 
+            Requires """
+        return
+
+    def test_nan_base(self):
+        return
+
+    def test_nan_values(self):
+        return
+
+    def blanding_nan_values(self):
+        return 
+
+    
+    def test_zero_halfwidth(self):
+        with self.assertRaises(ValueError) as e:
+            gridpp.calc_gradient(self.input_base, self.input_values, -1)
+    
+    
+    def test_negative_halfwidth(self):
+        return
+
+
+if __name__ == '__main__':
+	unittest.main()
+
+
+"""
+Notes
+
+self.input_base = [[0, 0, 0], [1, 1, 1], [4, 4, 4]]
+self.input_base[0:2, 0:2] = np.nan
+
+input size is null, return empty vector
+
+check test_invalid arguments in neighbourhood_test.py line 31
+
+use test with squared temperatures 
+"""
\ No newline at end of file