From caa19a9f456bd13d96f126c44f4c9279acb9817c Mon Sep 17 00:00:00 2001
From: Anton Afanasyev <anstaf@gmail.com>
Date: Tue, 28 Sep 2021 11:56:21 +0200
Subject: [PATCH 1/3] skip unit stride consistency check if the dimension size
 is 1

---
 .../storage/adapter/python_sid_adapter.hpp    | 36 +++++++++++++------
 1 file changed, 26 insertions(+), 10 deletions(-)

diff --git a/include/gridtools/storage/adapter/python_sid_adapter.hpp b/include/gridtools/storage/adapter/python_sid_adapter.hpp
index 93c382bc4..faa1c197d 100644
--- a/include/gridtools/storage/adapter/python_sid_adapter.hpp
+++ b/include/gridtools/storage/adapter/python_sid_adapter.hpp
@@ -22,7 +22,6 @@
 
 #include <pybind11/pybind11.h>
 
-#include <boost/optional.hpp>
 #include <boost/variant.hpp>
 
 #include "../../common/array.hpp"
@@ -39,11 +38,12 @@ namespace gridtools {
         template <size_t, class>
         struct kind {};
 
-        template <int UnitStrideDim>
+        template <size_t UnitStrideDim>
         struct transform_strides_f {
+            bool m_unit_stride_can_be_used;
             template <size_t I, class T, std::enable_if_t<I == UnitStrideDim, int> = 0>
             integral_constant<pybind11::ssize_t, 1> operator()(T val) const {
-                if (val != 1)
+                if (m_unit_stride_can_be_used && val != 1)
                     throw std::domain_error("incompatible strides, expected unit stride");
                 return {};
             }
@@ -54,7 +54,25 @@ namespace gridtools {
             }
         };
 
-        template <class T, size_t Dim, class Kind, int UnitStrideDim>
+        template <size_t UnitStrideDim,
+            class Strides,
+            class Shape,
+            std::enable_if_t<(UnitStrideDim >= tuple_util::size<std::decay_t<Strides>>::value), int> = 0>
+        Strides &&assign_unit_stride(Strides &&strides, Shape &&) {
+            return std::forward<Strides>(strides);
+        }
+
+        template <size_t UnitStrideDim,
+            class Strides,
+            class Shape,
+            std::enable_if_t<(UnitStrideDim < tuple_util::size<std::decay_t<Strides>>::value), int> = 0>
+        decltype(auto) assign_unit_stride(Strides &&strides, Shape &&shape) {
+            return tuple_util::transform_index(
+                transform_strides_f<UnitStrideDim>{tuple_util::get<UnitStrideDim>(shape) > 1},
+                tuple_util::convert_to<tuple>(std::forward<Strides>(strides)));
+        }
+
+        template <class T, size_t Dim, class Kind, size_t UnitStrideDim>
         struct wrapper {
             pybind11::buffer_info m_info;
 
@@ -68,8 +86,7 @@ namespace gridtools {
                     assert(obj.m_info.strides[i] % obj.m_info.itemsize == 0);
                     res[i] = obj.m_info.strides[i] / obj.m_info.itemsize;
                 }
-                return tuple_util::transform_index(
-                    transform_strides_f<UnitStrideDim>{}, tuple_util::convert_to<tuple>(res));
+                return assign_unit_stride<UnitStrideDim>(std::move(res), sid_get_upper_bounds(obj));
             }
 
             friend std::array<integral_constant<pybind11::ssize_t, 0>, Dim> sid_get_lower_bounds(wrapper const &) {
@@ -90,7 +107,7 @@ namespace gridtools {
             }
         };
 
-        template <class T, std::size_t Dim, class Kind = void, int UnitStrideDim = -1>
+        template <class T, std::size_t Dim, class Kind = void, size_t UnitStrideDim = size_t(-1)>
         wrapper<T, Dim, Kind, UnitStrideDim> as_sid(pybind11::buffer const &src) {
             static_assert(
                 std::is_trivially_copyable<T>::value, "as_sid should be instantiated with the trivially copyable type");
@@ -243,7 +260,7 @@ namespace gridtools {
             return res;
         }
 
-        template <class T, size_t Dim, class Kind = void, int UnitStrideDim = -1>
+        template <class T, size_t Dim, class Kind = void, size_t UnitStrideDim = size_t(-1)>
         auto as_cuda_sid(pybind11::object const &src) {
             static_assert(std::is_trivially_copyable<T>::value,
                 "as_cuda_sid should be instantiated with the trivially copyable type");
@@ -320,8 +337,7 @@ namespace gridtools {
             using sid::property;
             return sid::synthetic()
                 .template set<property::origin>(sid::host_device::simple_ptr_holder<T *>{ptr})
-                .template set<property::strides>(tuple_util::transform_index(
-                    transform_strides_f<UnitStrideDim>{}, tuple_util::convert_to<tuple>(strides)))
+                .template set<property::strides>(assign_unit_stride<UnitStrideDim>(std::move(strides), shape))
                 .template set<property::strides_kind, kind<Dim, Kind>>()
                 .template set<property::lower_bounds>(array<integral_constant<size_t, 0>, Dim>())
                 .template set<property::upper_bounds>(shape);

From bb98d7d5cb9d1612f95b112f491bf20fc9727d62 Mon Sep 17 00:00:00 2001
From: Anton Afanasyev <anstaf@gmail.com>
Date: Tue, 28 Sep 2021 12:03:44 +0200
Subject: [PATCH 2/3] add missing std::forward

---
 include/gridtools/storage/adapter/python_sid_adapter.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/gridtools/storage/adapter/python_sid_adapter.hpp b/include/gridtools/storage/adapter/python_sid_adapter.hpp
index faa1c197d..11b0b8f1d 100644
--- a/include/gridtools/storage/adapter/python_sid_adapter.hpp
+++ b/include/gridtools/storage/adapter/python_sid_adapter.hpp
@@ -68,7 +68,7 @@ namespace gridtools {
             std::enable_if_t<(UnitStrideDim < tuple_util::size<std::decay_t<Strides>>::value), int> = 0>
         decltype(auto) assign_unit_stride(Strides &&strides, Shape &&shape) {
             return tuple_util::transform_index(
-                transform_strides_f<UnitStrideDim>{tuple_util::get<UnitStrideDim>(shape) > 1},
+                transform_strides_f<UnitStrideDim>{tuple_util::get<UnitStrideDim>(std::forward<Shape>(shape)) > 1},
                 tuple_util::convert_to<tuple>(std::forward<Strides>(strides)));
         }
 

From ae698138bc7171e90f9610a258581aa6df2244f6 Mon Sep 17 00:00:00 2001
From: Anton Afanasyev <anstaf@gmail.com>
Date: Wed, 29 Sep 2021 10:38:47 +0200
Subject: [PATCH 3/3] add comment

---
 include/gridtools/storage/adapter/python_sid_adapter.hpp | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/include/gridtools/storage/adapter/python_sid_adapter.hpp b/include/gridtools/storage/adapter/python_sid_adapter.hpp
index 11b0b8f1d..245d5ce3c 100644
--- a/include/gridtools/storage/adapter/python_sid_adapter.hpp
+++ b/include/gridtools/storage/adapter/python_sid_adapter.hpp
@@ -67,8 +67,12 @@ namespace gridtools {
             class Shape,
             std::enable_if_t<(UnitStrideDim < tuple_util::size<std::decay_t<Strides>>::value), int> = 0>
         decltype(auto) assign_unit_stride(Strides &&strides, Shape &&shape) {
-            return tuple_util::transform_index(
-                transform_strides_f<UnitStrideDim>{tuple_util::get<UnitStrideDim>(std::forward<Shape>(shape)) > 1},
+            // Numpy may shuffle array layout if the array sizes are equal to one in some dimensions.
+            // In this case the static unit stride calculation doesn't match the actual strides.
+            // Luckily we can ignore that because those strides will not be used (corresponding shape == 1).
+            // See: https://numpy.org/devdocs/release/1.8.0-notes.html#npy-relaxed-strides-checking
+            bool unit_stride_can_be_used = tuple_util::get<UnitStrideDim>(std::forward<Shape>(shape)) > 1;
+            return tuple_util::transform_index(transform_strides_f<UnitStrideDim>{unit_stride_can_be_used},
                 tuple_util::convert_to<tuple>(std::forward<Strides>(strides)));
         }