Repository: xtensor-stack/xtensor (language: python)

Issue:
bug in histogram function
Hi, I'm noticing unexpected behavior in `xt::histogram()`. In particular, given some `data` that exceeds the boundaries in `bin_edges`, the data is put into the last bin. I believe this is because in `xhistogram.hpp`, the range is set incorrectly here: `while (data[idx] >= bin_edges[ibin + 1] && ibin < bin_edges.size() - 2)`. When bin_edges has size four, for instance, the loop doesn't run after `ibin` is greater than or equal to `bin_edges.size() - 2 = 2`, and hence, the last bin (bin 2 in this case) accumulates the rest of the data: `count[ibin] += weights[idx];`. Relevant screenshots are attached:
<img width="1140" alt="Screen Shot 2020-07-23 at 4 33 25 PM" src="https://user-images.githubusercontent.com/29467219/88337038-04a97880-cd04-11ea-9686-89f1fd21dfc9.png">
<img width="1140" alt="Screen Shot 2020-07-23 at 4 33 43 PM" src="https://user-images.githubusercontent.com/29467219/88337054-0d01b380-cd04-11ea-8f95-ec7c136f4465.png">

This can be resolved by instead using terminal condition `while (data[idx] >= bin_edges[ibin + 1] && ibin < bin_edges.size() - 1)`, and dropping the last bin. Not sure if this is the anticipated behavior, but it deviates from the NumPy result.


Discussion / hints:
Thanks for reporting.  This is indeed not the expected behavior. Notice there is [oin-going work](https://github.com/xtensor-stack/xtensor/pull/2088) on the histogram function, so maybe this could also be fixed in this PR?
cc @zhujun98
I am back! Sorry for having disappeared for a long time. Yes, this can be solved in #2088 for histogram with equally sized bins. But another PR will need to address the issue for the rest `histogram_algorithm`s.

Fix patch:
diff --git a/include/xtensor/xhistogram.hpp b/include/xtensor/xhistogram.hpp
index 08c1ea9b8..7b8f75d4f 100644
--- a/include/xtensor/xhistogram.hpp
+++ b/include/xtensor/xhistogram.hpp
@@ -58,8 +58,6 @@ namespace xt
             XTENSOR_ASSERT(weights.size() == data.size());
             XTENSOR_ASSERT(bin_edges.size() >= 2);
             XTENSOR_ASSERT(std::is_sorted(bin_edges.cbegin(), bin_edges.cend()));
-            XTENSOR_ASSERT(xt::amin(data)[0] >= bin_edges[0]);
-            XTENSOR_ASSERT(xt::amax(data)[0] <= bin_edges[bin_edges.size() - 1]);
 
             size_t n_bins = bin_edges.size() - 1;
             xt::xtensor<value_type, 1> count = xt::zeros<value_type>({ n_bins });
@@ -93,10 +91,23 @@ namespace xt
 
                 for (auto& idx : sorter)
                 {
-                    while (data[idx] >= bin_edges[ibin + 1] && ibin < bin_edges.size() - 2)
+                    auto item = data[idx];
+
+                    if (item < bin_edges[0])
+                    {
+                        continue;
+                    }
+
+                    if (item > bin_edges[n_bins])
+                    {
+                        break;
+                    }
+
+                    while (item >= bin_edges[ibin + 1] && ibin < n_bins - 1)
                     {
                         ++ibin;
                     }
+
                     count[ibin] += weights[idx];
                 }
             }
@@ -247,8 +258,7 @@ namespace xt
         // - size
         XTENSOR_ASSERT(weights.size() == data.size());
         // - bounds
-        XTENSOR_ASSERT(left <= xt::amin(data)[0]);
-        XTENSOR_ASSERT(right >= xt::amax(data)[0]);
+        XTENSOR_ASSERT(left <= right);
         // - non-empty
         XTENSOR_ASSERT(bins > std::size_t(0));
 


Test patch:
diff --git a/test/files/cppy_source/test_extended_broadcast_view.cppy b/test/files/cppy_source/test_extended_broadcast_view.cppy
index 5aab7afb0..0852d5e2a 100644
--- a/test/files/cppy_source/test_extended_broadcast_view.cppy
+++ b/test/files/cppy_source/test_extended_broadcast_view.cppy
@@ -1,5 +1,6 @@
 /***************************************************************************
-* Copyright (c) 2016, Johan Mabille, Sylvain Corlay and Wolf Vollprecht    *
+* Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
+* Copyright (c) QuantStack                                                 *
 *                                                                          *
 * Distributed under the terms of the BSD 3-Clause License.                 *
 *                                                                          *
diff --git a/test/files/cppy_source/test_extended_xhistogram.cppy b/test/files/cppy_source/test_extended_xhistogram.cppy
index 5b32915d4..61139be9d 100644
--- a/test/files/cppy_source/test_extended_xhistogram.cppy
+++ b/test/files/cppy_source/test_extended_xhistogram.cppy
@@ -1,5 +1,6 @@
 /***************************************************************************
-* Copyright (c) 2016, Johan Mabille, Sylvain Corlay and Wolf Vollprecht    *
+* Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
+* Copyright (c) QuantStack                                                 *
 *                                                                          *
 * Distributed under the terms of the BSD 3-Clause License.                 *
 *                                                                          *
@@ -45,4 +46,26 @@ namespace xt
 
         EXPECT_EQ(prob, py_prob);
     }
+
+    /*py
+
+    a = (np.random.random(1000) * 51).astype(np.int)
+
+    count, bin_edges  = np.histogram(a, 30, (10, 40))
+
+    */
+    TEST(xtest_extended_histogram, histogram_with_range)
+    {
+        // py_a
+        // py_bin_edges
+        // py_count
+
+        auto bin_edges = xt::histogram_bin_edges(py_a, 10, 40, std::size_t(30));
+
+        auto count = xt::histogram(py_a, bin_edges);
+
+        EXPECT_EQ(xt::cast<double>(bin_edges), py_bin_edges);
+
+        EXPECT_EQ(count, xt::cast<double>(py_count));
+    }
 }
diff --git a/test/files/cppy_source/test_extended_xmath_interp.cppy b/test/files/cppy_source/test_extended_xmath_interp.cppy
index 6058e55d7..0dc17aa7a 100644
--- a/test/files/cppy_source/test_extended_xmath_interp.cppy
+++ b/test/files/cppy_source/test_extended_xmath_interp.cppy
@@ -1,5 +1,6 @@
 /***************************************************************************
-* Copyright (c) 2016, Johan Mabille, Sylvain Corlay and Wolf Vollprecht    *
+* Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
+* Copyright (c) QuantStack                                                 *
 *                                                                          *
 * Distributed under the terms of the BSD 3-Clause License.                 *
 *                                                                          *
diff --git a/test/files/cppy_source/test_extended_xmath_reducers.cppy b/test/files/cppy_source/test_extended_xmath_reducers.cppy
index 4160f1bee..13aedd3ee 100644
--- a/test/files/cppy_source/test_extended_xmath_reducers.cppy
+++ b/test/files/cppy_source/test_extended_xmath_reducers.cppy
@@ -1,3 +1,11 @@
+/***************************************************************************
+* Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
+* Copyright (c) QuantStack                                                 *
+*                                                                          *
+* Distributed under the terms of the BSD 3-Clause License.                 *
+*                                                                          *
+* The full license is in the file LICENSE, distributed with this software. *
+****************************************************************************/
 
 #include <algorithm>
 
diff --git a/test/files/cppy_source/test_extended_xsort.cppy b/test/files/cppy_source/test_extended_xsort.cppy
index f477a0dfe..187cc059f 100644
--- a/test/files/cppy_source/test_extended_xsort.cppy
+++ b/test/files/cppy_source/test_extended_xsort.cppy
@@ -1,5 +1,6 @@
 /***************************************************************************
-* Copyright (c) 2016, Johan Mabille, Sylvain Corlay and Wolf Vollprecht    *
+* Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
+* Copyright (c) QuantStack                                                 *
 *                                                                          *
 * Distributed under the terms of the BSD 3-Clause License.                 *
 *                                                                          *
diff --git a/test/files/preprocess.py b/test/files/preprocess.py
index b3a21abbd..adf02d22f 100644
--- a/test/files/preprocess.py
+++ b/test/files/preprocess.py
@@ -68,7 +68,10 @@ def exec_comment(txt, upper_level=False):
 	for line in contents.split('\n'):
 		if lidx == 8:
 			f = os.path.split(f)[1]
-			result_file += "// This file is generated from test/files/cppy_source/{} by preprocess.py!".format(f) + '\n\n'
+			result_file += "// This file is generated from the source file test/files/cppy_source/{} by preprocess.py!\n".format(f)
+			result_file += "// Warning: This file should not be modified directly! " \
+						   "Please modify the source file to change the test!\n\n"
+
 		lstrip = line.lstrip()
 		if lstrip.startswith("/*py"):
 			exec_comment(matches[idx], True)
diff --git a/test/test_extended_xhistogram.cpp b/test/test_extended_xhistogram.cpp
index 78d882bec..31f543620 100644
--- a/test/test_extended_xhistogram.cpp
+++ b/test/test_extended_xhistogram.cpp
@@ -6,8 +6,9 @@
 *                                                                          *
 * The full license is in the file LICENSE, distributed with this software. *
 ****************************************************************************/
+// This file is generated from the source file test/files/cppy_source/test_extended_xhistogram.cppy by preprocess.py!
+// Warning: This file should not be modified directly! Please modify the source file to change the test!
 
-// This file is generated from test/files/cppy_source/test_extended_xhistogram.cppy by preprocess.py!
 
 #include <algorithm>
 
@@ -101,5 +102,73 @@ namespace xt
 
         EXPECT_EQ(prob, py_prob);
     }
+
+    /*py
+
+    a = (np.random.random(1000) * 51).astype(np.int)
+
+    count, bin_edges  = np.histogram(a, 30, (10, 40))
+
+    */
+    TEST(xtest_extended_histogram, histogram_with_range)
+    {
+        // py_a
+        xarray<long> py_a = { 9,27,44,37,41,33,35,43,12,24,11,50,48, 2,35,47, 9,28,46, 1,35,15,47,49,
+                             48,24,43,43,16,42, 1,30,11, 6, 3,35,17,36, 3,16,27,40,16,31,45,31,11, 1,
+                             44, 1,44,26,47,40,50,17,39,20,24,32,44,50,39,21,21,37,12, 5,18,14,15,11,
+                              2, 0,50,21,19,34,11,48,40, 4,21,44,48,23,31, 8,50,11,48,33,30,26,11, 9,
+                             11, 9,39,17, 2,49,45,47,50, 8,20,38,35, 7,41,11,11,27,30,29, 4,44,13, 6,
+                             45,48,43,41,33,28, 4,20,19,13,36,25, 4,11,34, 3,43,25,24,30,42,17,34,28,
+                             13,44,40,33,43,44,36,42,35,34,31,38, 8,44,44, 1,42, 6,17,37, 8,41,42,25,
+                              0,14,31,50,32,13,32,27,39, 5,38,27,49,17,32,47, 5,47,35, 3,15,36, 3,29,
+                             17,31, 2,44,49,49,38, 6,38, 1, 1,16,24,39,34,22,13,50,21,23, 8,40,35,11,
+                              4,34,33,13,48, 7,22,48,21,32,20,13,50,20,45,11,10, 1,33,18,44,24,49, 9,
+                             44,39,39,43,38,31, 6, 1,46,31,40,24, 5, 6,34,21,10,25, 3,29,13,40,15,23,
+                              0, 3,20,24,30,14,35,43,39, 2,24, 5,12,50, 7,25,31,35,28, 0,16,26, 4,17,
+                              1, 4,20, 6,28,35,40,10, 8, 5,32,36, 1,47, 2,27,36,44,36,40,17,41, 4,45,
+                             27,41,23,32,26,37, 4, 3,12, 8,44,11,49,17, 9,40,33,25,28,36,11,50,49,33,
+                             10,34, 3, 1,13,23,44,37,37,21,17,18,50, 2,44,29,22,36,24,44,45,21,14,30,
+                             46,10,31,32,37, 6,36,46, 9,12,49, 9,43,25,12,44,22,26,18,30, 8,19,49,13,
+                             33,16,39, 6,49,23,12, 3, 8,26,17,42,21,12,31,36, 8, 8, 1,37,33,24,43,41,
+                             29,44,10, 5,13, 2,27,47, 2, 6,23,47,16,25, 2, 7,50,49, 0,48,32,44,23,26,
+                             24,34, 7, 1,15,35,10,34,49, 4,34,22,44, 9,35,42,48,34,25,31,44,29, 1,47,
+                             35,34,10,33,20,33, 5,33,50, 2,49,20,44,39,28,37,44,20,16,34,41,38,40,22,
+                             41, 6,27, 0,16,18,20,35,19,22,12,19,11, 3,30,34,31,23,19,44,26,24, 1,17,
+                             19,20,29,27,31,39,41,36,48, 0, 9, 0,33,45,12,47, 3,47,17, 5,24,13,14,15,
+                             40,27,15,31,36,13,21, 6, 9,34, 9,26,36, 5,28,13,49,24,41,28, 2,32,48,30,
+                             41,45,11,10,31,20,42,45,18,12,39,14,41,21,34, 4,31,23,29, 8,37,44,11, 4,
+                              1,32,30,27,11,19,30,25,50, 6,35,20,21,36,35,50, 6, 5,36,29,13, 4, 4,45,
+                              9,16,11,18, 3,26, 3,40,11,27,44,33,27,16,16,34,50,33,28,37,23, 3,28,48,
+                              8,35,10,27, 4,22,38,17,33,40,47,11,20, 7,50,47,27,42,26,31, 4,38, 6,42,
+                             39,36, 1,15,13,18, 4,47,28,15,20,22,30,26,46,25,50,43,10,47, 5,41,19,44,
+                             44,41,40,15, 4,20, 8,35,17,49,32,41, 6,43,47,24,30,39, 8,25,20, 7,18, 3,
+                              1, 6,49,28,49,22,15,25,22, 5,32,11,31,33, 7, 3,39,23, 2,50, 2,35,50,12,
+                              7, 6,15, 5,35, 3,25,50,41,31,15,31,26,21, 6,45,22, 9,18,21,42,37,39, 0,
+                             21,24, 0,13,38, 6,27,10, 0,12,49,40,48,24, 5,27,23,43, 5,24, 7,16,37,24,
+                             19,20,23,40,45,48,40,16,35,22,12,42, 1,45,23,32,33,45,32,31, 3,26, 7,37,
+                             26,34, 2, 4,36, 3, 3, 0,48,37,18,15,17,39,33, 9, 8, 5,33,38,13, 1, 4,49,
+                             15,39,31,19,10, 6,31,39,32,27, 2,49,40,14,49,30,29,38,41,33, 6,17,47,11,
+                             18,22,22,31,48,12, 6,10,45,32,14,41,43,43,46,12,38,23,42,37,39,33, 9,27,
+                             50,47, 2, 8, 6,37,41,10,25,42,37,27,30,25,15,28,35,44,32,38, 8,23, 0,12,
+                             37,50, 5,20,40,10,28,37,31, 9,18,39,28, 0,38, 1,38,10,48,18,16, 7,15,44,
+                             50,18,22,36,45,30,19,21,35, 0,31,18,40, 4,29,24,32, 3,29,28,28,30,34,41,
+                             13,42,25, 3, 2,17,40,36,40,26,22, 7,16,22, 4,11,30,37,50,47,32,21,32,40,
+                              6,20,42,19,29,29, 9,18,17, 1, 1,42,13,26,15,47,13,21,44,42, 9,40,23,24,
+                              6, 4,37,25,22,37,39, 8,31, 6,38,33,48, 3, 2,14};
+        // py_bin_edges
+        xarray<double> py_bin_edges = {10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,
+                                       28.,29.,30.,31.,32.,33.,34.,35.,36.,37.,38.,39.,40.};
+        // py_count
+        xarray<long> py_count = {17,25,17,21, 9,18,16,20,17,13,21,19,19,18,22,18,17,21,18,14,18,28,21,24,
+                                 20,24,20,23,17,47};
+
+        auto bin_edges = xt::histogram_bin_edges(py_a, 10, 40, std::size_t(30));
+
+        auto count = xt::histogram(py_a, bin_edges);
+
+        EXPECT_EQ(xt::cast<double>(bin_edges), py_bin_edges);
+
+        EXPECT_EQ(count, xt::cast<double>(py_count));
+    }
 }