src/cpu/kernels/pool3d/neon/impl.cpp - ml/ComputeLibrary - Gitiles

 /*
  * Copyright (c) 2022 Arm Limited.
  *
  * SPDX-License-Identifier: MIT
  *
  * Permission is hereby granted, free of charge, to any person obtaining a copy
  * of this software and associated documentation files (the "Software"), to
  * deal in the Software without restriction, including without limitation the
  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
  * sell copies of the Software, and to permit persons to whom the Software is
  * furnished to do so, subject to the following conditions:
  *
  * The above copyright notice and this permission notice shall be included in all
  * copies or substantial portions of the Software.
  *
  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  * SOFTWARE.
  */
 #include "arm_compute/core/Helpers.h"
 #include "arm_compute/core/ITensor.h"
 #include "arm_compute/core/Types.h"
 #include "arm_compute/core/utils/misc/Traits.h"
 #include "src/core/NEON/wrapper/intrinsics/intrinsics.h"
 #include "src/core/helpers/WindowHelpers.h"

 #include "src/cpu/kernels/pool3d/neon/impl.h"

 namespace arm_compute
 {
 namespace cpu
 {
 namespace
 {
 inline float calculate_avg_scale(bool exclude_padding, const Coordinates &id, const int pool_size_x, const int pool_size_y, const int pool_size_z, const int upper_bound_w,
                                  const int upper_bound_h, const int upper_bound_d, const int pad_x, const int pad_y, const int pad_z, const int stride_x, const int stride_y, const int stride_z)
 {
     // Based on NDHWC
     int start_x = id[1] * stride_x - pad_x;
     int start_y = id[2] * stride_y - pad_y;
     int start_z = id[3] * stride_z - pad_z;

     const int end_x = std::min(start_x + pool_size_x, upper_bound_w);
     const int end_y = std::min(start_y + pool_size_y, upper_bound_h);
     const int end_z = std::min(start_z + pool_size_z, upper_bound_d);
     if(exclude_padding)
     {
         start_x = std::max(0, start_x);
         start_y = std::max(0, start_y);
         start_z = std::max(0, start_z);
     }
     return 1.f / ((end_y - start_y) * (end_x - start_x) * (end_z - start_z));
 }


 template <typename T>
 void max_poolingMxNxD_fp_neon_ndhwc(const ITensor *src, ITensor *dst0, Pooling3dLayerInfo &pool_info, const Window &window_out,
                                     const int window_start_x, const int window_end_x, const int window_step_x)

 {
     using vtype       = wrapper::traits::neon_bitvector<T, wrapper::traits::BitWidth::W128>;
     using vector_type = typename vtype::type;
     using tag_type    = typename vtype::tag_type;

     int pool_stride_x = static_cast<int>(pool_info.stride.width);
     int pool_stride_y = static_cast<int>(pool_info.stride.height);
     int pool_stride_z = static_cast<int>(pool_info.stride.depth);

     const int pool_size_x = pool_info.is_global_pooling ? src->info()->tensor_shape().y() : pool_info.pool_size.width;
     const int pool_size_y = pool_info.is_global_pooling ? src->info()->tensor_shape().z() : pool_info.pool_size.height;
     const int pool_size_z = pool_info.is_global_pooling ? src->info()->tensor_shape()[3] : pool_info.pool_size.depth;

     const int pool_pad_top   = static_cast<int>(pool_info.padding.top);
     const int pool_pad_left  = static_cast<int>(pool_info.padding.left);
     const int pool_pad_front = static_cast<int>(pool_info.padding.front);

     const int input_dim_w = src->info()->dimension(1);
     const int input_dim_h = src->info()->dimension(2);
     const int input_dim_d = src->info()->dimension(3);

     const int y_stride = static_cast<int>(src->info()->strides_in_bytes().y());
     const int z_stride = static_cast<int>(src->info()->strides_in_bytes().z());
     const int w_stride = static_cast<int>(src->info()->strides_in_bytes()[3]);
     const int n_stride = static_cast<int>(src->info()->strides_in_bytes()[4]);

     const uint8_t *in_ptr_start = src->buffer() + src->info()->offset_first_element_in_bytes();

     Iterator out(dst0, window_out);

     vector_type vres;
     execute_window_loop(window_out, [&](const Coordinates & id)
     {
         // Computing the theoretical input starting/ending points
         const int in_idx_width  = static_cast<int>(id.y()) * pool_stride_x - pool_pad_left;
         const int in_idx_height = static_cast<int>(id.z()) * pool_stride_y - pool_pad_top;
         const int in_idx_depth  = static_cast<int>(id[3]) * pool_stride_z - pool_pad_front;

         const int pool_start_x = std::max(0, -in_idx_width);
         const int pool_end_x_t = std::min(input_dim_w + pool_pad_left - in_idx_width, pool_size_x);
         const int pool_start_y = std::max(0, -in_idx_height);
         const int pool_end_y_t = std::min(input_dim_h + pool_pad_top - in_idx_height, pool_size_y);

         const int pool_start_z = std::max(0, -in_idx_depth);
         const int pool_end_z_t = std::min(input_dim_d + pool_pad_front - in_idx_depth, pool_size_z);

         // The end of width to consider in calculation should exclude PAD_X, PAD_Y and PAD_Z
         const int pool_end_x = std::min(pool_end_x_t, input_dim_w - in_idx_width);
         const int pool_end_y = std::min(pool_end_y_t, input_dim_h - in_idx_height);
         const int pool_end_z = std::min(pool_end_z_t, input_dim_d - in_idx_depth);

         const uint8_t *in_ptr_n = in_ptr_start + id[4] * n_stride;

         int x_off = window_start_x;

         for(; x_off <= (window_end_x - window_step_x); x_off += window_step_x) // C
         {
             vres = wrapper::vdup_n(static_cast<T>(-std::numeric_limits<float>::infinity()), tag_type());
             for(int z = pool_start_z; z < pool_end_z; ++z)
             {
                 const uint8_t *in_ptr_z = in_ptr_n + (z + in_idx_depth) * w_stride;
                 for(int y = pool_start_y; y < pool_end_y; ++y)
                 {
                     const uint8_t *in_ptr_y = in_ptr_z + (y + in_idx_height) * z_stride;
                     for(int x = pool_start_x; x < pool_end_x; ++x)
                     {
                         const uint8_t    *in_ptr_x = in_ptr_y + (x + in_idx_width) * y_stride;
                         const vector_type data     = wrapper::vloadq(reinterpret_cast<const T *>(in_ptr_x) + x_off);
                         vres                       = wrapper::vmax(vres, data);
                     }
                 }
             }
             // Store result
             wrapper::vstore(reinterpret_cast<T *>(out.ptr()) + x_off, vres);
         }

         // Left-overs loop
         for(; x_off < window_end_x; ++x_off)
         {
             T res(0);
             res = -std::numeric_limits<float>::infinity();
             for(int z = pool_start_z; z < pool_end_z; ++z)
             {
                 const uint8_t *in_ptr_z = in_ptr_n + (z + in_idx_depth) * w_stride;
                 for(int y = pool_start_y; y < pool_end_y; ++y)
                 {
                     const uint8_t *in_ptr_y = in_ptr_z + (y + in_idx_height) * z_stride;
                     for(int x = pool_start_x; x < pool_end_x; ++x)
                     {
                         const uint8_t *in_ptr_x = in_ptr_y + (x + in_idx_width) * y_stride;
                         const T        data     = *(reinterpret_cast<const T *>(in_ptr_x) + x_off);
                         res                     = std::max(res, data);
                     }
                 }
             }
             // Store result
             *(reinterpret_cast<T *>(out.ptr()) + x_off) = res;
         }
     },
     out);
 }

 template <typename T>
 void avg_poolingMxNxD_fp_neon_ndhwc(const ITensor *src, ITensor *dst0, Pooling3dLayerInfo &pool_info,
                                     const Window &window_out, const int window_start_x, const int window_end_x, const int window_step_x)
 {
     using vtype       = wrapper::traits::neon_bitvector<T, wrapper::traits::BitWidth::W128>;
     using vector_type = typename vtype::type;
     using tag_type    = typename vtype::tag_type;

     int pool_stride_x = static_cast<int>(pool_info.stride.width);
     int pool_stride_y = static_cast<int>(pool_info.stride.height);
     int pool_stride_z = static_cast<int>(pool_info.stride.depth);

     const int pool_size_x = pool_info.is_global_pooling ? src->info()->tensor_shape().y() : pool_info.pool_size.width;
     const int pool_size_y = pool_info.is_global_pooling ? src->info()->tensor_shape().z() : pool_info.pool_size.height;
     const int pool_size_z = pool_info.is_global_pooling ? src->info()->tensor_shape()[3] : pool_info.pool_size.depth;

     const int pool_pad_top    = static_cast<int>(pool_info.padding.top);
     const int pool_pad_bottom = static_cast<int>(pool_info.padding.bottom);
     const int pool_pad_left   = static_cast<int>(pool_info.padding.left);
     const int pool_pad_right  = static_cast<int>(pool_info.padding.right);
     const int pool_pad_front  = static_cast<int>(pool_info.padding.front);
     const int pool_pad_back   = static_cast<int>(pool_info.padding.back);

     const int upper_bound_w = src->info()->dimension(1) + (pool_info.exclude_padding ? 0 : pool_pad_right);
     const int upper_bound_h = src->info()->dimension(2) + (pool_info.exclude_padding ? 0 : pool_pad_bottom);
     const int upper_bound_d = src->info()->dimension(3) + (pool_info.exclude_padding ? 0 : pool_pad_back);

     const int input_dim_w = src->info()->dimension(1);
     const int input_dim_h = src->info()->dimension(2);
     const int input_dim_d = src->info()->dimension(3);

     const int y_stride = static_cast<int>(src->info()->strides_in_bytes().y());
     const int z_stride = static_cast<int>(src->info()->strides_in_bytes().z());
     const int w_stride = static_cast<int>(src->info()->strides_in_bytes()[3]);
     const int n_stride = static_cast<int>(src->info()->strides_in_bytes()[4]);

     const uint8_t *in_ptr_start = src->buffer() + src->info()->offset_first_element_in_bytes();

     Iterator out(dst0, window_out);

     vector_type vres;
     execute_window_loop(window_out, [&](const Coordinates & id)
     {
         // Computing the theoretical input starting/ending points
         const int in_idx_width  = static_cast<int>(id.y()) * pool_stride_x - pool_pad_left;
         const int in_idx_height = static_cast<int>(id.z()) * pool_stride_y - pool_pad_top;
         const int in_idx_depth  = static_cast<int>(id[3]) * pool_stride_z - pool_pad_front;

         const int pool_start_x = std::max(0, -in_idx_width);
         const int pool_end_x_t = std::min(input_dim_w + pool_pad_left - in_idx_width, pool_size_x);
         const int pool_start_y = std::max(0, -in_idx_height);
         const int pool_end_y_t = std::min(input_dim_h + pool_pad_top - in_idx_height, pool_size_y);

         const int pool_start_z = std::max(0, -in_idx_depth);
         const int pool_end_z_t = std::min(input_dim_d + pool_pad_front - in_idx_depth, pool_size_z);

         // The end of width to consider in calculation should exclude PAD_X, PAD_Y and PAD_Z
         const int pool_end_x = std::min(pool_end_x_t, input_dim_w - in_idx_width);
         const int pool_end_y = std::min(pool_end_y_t, input_dim_h - in_idx_height);
         const int pool_end_z = std::min(pool_end_z_t, input_dim_d - in_idx_depth);

         const uint8_t *in_ptr_n = in_ptr_start + id[4] * n_stride;

         // Calculate scale
         const float scale = calculate_avg_scale(pool_info.exclude_padding, id, pool_size_x, pool_size_y, pool_size_z, upper_bound_w, upper_bound_h, upper_bound_d, pool_pad_left,
                                                 pool_pad_top, pool_pad_front, pool_stride_x,
                                                 pool_stride_y, pool_stride_z);
         const vector_type scale_v = wrapper::vdup_n(static_cast<T>(scale), tag_type());

         int x_off = window_start_x;

         for(; x_off <= (window_end_x - window_step_x); x_off += window_step_x) // C
         {
             // Perform pooling
             vres = wrapper::vdup_n(static_cast<T>(0.0f), tag_type());
             for(int z = pool_start_z; z < pool_end_z; ++z)
             {
                 const uint8_t *in_ptr_z = in_ptr_n + (z + in_idx_depth) * w_stride;
                 for(int y = pool_start_y; y < pool_end_y; ++y)
                 {
                     const uint8_t *in_ptr_y = in_ptr_z + (y + in_idx_height) * z_stride;
                     for(int x = pool_start_x; x < pool_end_x; ++x)
                     {
                         const uint8_t    *in_ptr_x = in_ptr_y + (x + in_idx_width) * y_stride;
                         const vector_type data     = wrapper::vloadq(reinterpret_cast<const T *>(in_ptr_x) + x_off);
                         vres                       = wrapper::vadd(vres, data);
                     }
                 }
             }

             // Divide by scale
             vres = wrapper::vmul(vres, scale_v);

             // Store result
             wrapper::vstore(reinterpret_cast<T *>(out.ptr()) + x_off, vres);
         }

         // Left-overs loop
         for(; x_off < window_end_x; ++x_off)
         {
             T res(0);

             for(int z = pool_start_z; z < pool_end_z; ++z)
             {
                 const uint8_t *in_ptr_z = in_ptr_n + (z + in_idx_depth) * w_stride;
                 for(int y = pool_start_y; y < pool_end_y; ++y)
                 {
                     const uint8_t *in_ptr_y = in_ptr_z + (y + in_idx_height) * z_stride;
                     for(int x = pool_start_x; x < pool_end_x; ++x)
                     {
                         const uint8_t *in_ptr_x = in_ptr_y + (x + in_idx_width) * y_stride;
                         const T        data     = *(reinterpret_cast<const T *>(in_ptr_x) + x_off);
                         res += data;
                     }
                 }
             }

             // Divide by scale
             res *= scale;

             // Store result
             *(reinterpret_cast<T *>(out.ptr()) + x_off) = res;
         }
     },
     out);
 }

 template <typename T>
 void l2_poolingMxNxD_fp_neon_ndhwc(const ITensor *src, ITensor *dst0, Pooling3dLayerInfo &pool_info,
                                    const Window &window_out, const int window_start_x, const int window_end_x, const int window_step_x)
 {
     using vtype       = wrapper::traits::neon_bitvector<T, wrapper::traits::BitWidth::W128>;
     using vector_type = typename vtype::type;
     using tag_type    = typename vtype::tag_type;

     int pool_stride_x = static_cast<int>(pool_info.stride.width);
     int pool_stride_y = static_cast<int>(pool_info.stride.height);
     int pool_stride_z = static_cast<int>(pool_info.stride.depth);

     const int pool_size_x = pool_info.is_global_pooling ? src->info()->tensor_shape().y() : pool_info.pool_size.width;
     const int pool_size_y = pool_info.is_global_pooling ? src->info()->tensor_shape().z() : pool_info.pool_size.height;
     const int pool_size_z = pool_info.is_global_pooling ? src->info()->tensor_shape()[3] : pool_info.pool_size.depth;

     const int pool_pad_top    = static_cast<int>(pool_info.padding.top);
     const int pool_pad_bottom = static_cast<int>(pool_info.padding.bottom);
     const int pool_pad_left   = static_cast<int>(pool_info.padding.left);
     const int pool_pad_right  = static_cast<int>(pool_info.padding.right);
     const int pool_pad_front  = static_cast<int>(pool_info.padding.front);
     const int pool_pad_back   = static_cast<int>(pool_info.padding.back);

     const int upper_bound_w = src->info()->dimension(1) + (pool_info.exclude_padding ? 0 : pool_pad_right);
     const int upper_bound_h = src->info()->dimension(2) + (pool_info.exclude_padding ? 0 : pool_pad_bottom);
     const int upper_bound_d = src->info()->dimension(3) + (pool_info.exclude_padding ? 0 : pool_pad_back);

     const int input_dim_w = src->info()->dimension(1);
     const int input_dim_h = src->info()->dimension(2);
     const int input_dim_d = src->info()->dimension(3);

     const int y_stride = static_cast<int>(src->info()->strides_in_bytes().y());
     const int z_stride = static_cast<int>(src->info()->strides_in_bytes().z());
     const int w_stride = static_cast<int>(src->info()->strides_in_bytes()[3]);
     const int n_stride = static_cast<int>(src->info()->strides_in_bytes()[4]);

     const uint8_t *in_ptr_start = src->buffer() + src->info()->offset_first_element_in_bytes();

     Iterator out(dst0, window_out);

     vector_type vres;
     execute_window_loop(window_out, [&](const Coordinates & id)
     {
         // Computing the theoretical input starting/ending points
         const int in_idx_width  = static_cast<int>(id.y()) * pool_stride_x - pool_pad_left;
         const int in_idx_height = static_cast<int>(id.z()) * pool_stride_y - pool_pad_top;
         const int in_idx_depth  = static_cast<int>(id[3]) * pool_stride_z - pool_pad_front;

         const int pool_start_x = std::max(0, -in_idx_width);
         const int pool_end_x_t = std::min(input_dim_w + pool_pad_left - in_idx_width, pool_size_x);
         const int pool_start_y = std::max(0, -in_idx_height);
         const int pool_end_y_t = std::min(input_dim_h + pool_pad_top - in_idx_height, pool_size_y);

         const int pool_start_z = std::max(0, -in_idx_depth);
         const int pool_end_z_t = std::min(input_dim_d + pool_pad_front - in_idx_depth, pool_size_z);

         // The end of width to consider in calculation should exclude PAD_X, PAD_Y and PAD_Z
         const int pool_end_x = std::min(pool_end_x_t, input_dim_w - in_idx_width);
         const int pool_end_y = std::min(pool_end_y_t, input_dim_h - in_idx_height);
         const int pool_end_z = std::min(pool_end_z_t, input_dim_d - in_idx_depth);

         const uint8_t *in_ptr_n = in_ptr_start + id[4] * n_stride;

         // Calculate scale
         const float scale = calculate_avg_scale(pool_info.exclude_padding, id, pool_size_x, pool_size_y, pool_size_z, upper_bound_w, upper_bound_h, upper_bound_d, pool_pad_left,
                                                 pool_pad_top, pool_pad_front, pool_stride_x,
                                                 pool_stride_y, pool_stride_z);

         int x_off = window_start_x;

         for(; x_off <= (window_end_x - window_step_x); x_off += window_step_x) // C
         {
             // Perform pooling
             vres = wrapper::vdup_n(static_cast<T>(0.0f), tag_type());
             for(int z = pool_start_z; z < pool_end_z; ++z)
             {
                 const uint8_t *in_ptr_z = in_ptr_n + (z + in_idx_depth) * w_stride;
                 for(int y = pool_start_y; y < pool_end_y; ++y)
                 {
                     const uint8_t *in_ptr_y = in_ptr_z + (y + in_idx_height) * z_stride;
                     for(int x = pool_start_x; x < pool_end_x; ++x)
                     {
                         const uint8_t    *in_ptr_x = in_ptr_y + (x + in_idx_width) * y_stride;
                         const vector_type data     = wrapper::vloadq(reinterpret_cast<const T *>(in_ptr_x) + x_off);
                         vres                       = wrapper::vmla(vres, data, data);
                     }
                 }
             }

             const vector_type scale_v = wrapper::vdup_n(static_cast<T>(scale), tag_type());

             // Divide by scale
             vres = wrapper::vmul(vres, scale_v);

             // Calculate square-root
             vres = wrapper::vinv(wrapper::vinvsqrt(vres));

             // Store result
             wrapper::vstore(reinterpret_cast<T *>(out.ptr()) + x_off, vres);
         }

         // Left-overs loop
         for(; x_off < window_end_x; ++x_off)
         {
             T res(0);

             for(int z = pool_start_z; z < pool_end_z; ++z)
             {
                 const uint8_t *in_ptr_z = in_ptr_n + (z + in_idx_depth) * w_stride;
                 for(int y = pool_start_y; y < pool_end_y; ++y)
                 {
                     const uint8_t *in_ptr_y = in_ptr_z + (y + in_idx_height) * z_stride;
                     for(int x = pool_start_x; x < pool_end_x; ++x)
                     {
                         const uint8_t *in_ptr_x = in_ptr_y + (x + in_idx_width) * y_stride;
                         const T        data     = *(reinterpret_cast<const T *>(in_ptr_x) + x_off);
                         res += data * data;
                     }
                 }
             }

             // Divide by scale
             res *= scale;

             // Square root
             res = std::sqrt(res);

             // Store result
             *(reinterpret_cast<T *>(out.ptr()) + x_off) = res;
         }
     },
     out);
 }
 } // namespace

 template <typename T>
 void poolingMxNxD_fp_neon_ndhwc(const ITensor *src, ITensor *dst0, Pooling3dLayerInfo &pool_info, const Window &window)
 {
     const int     window_start_x = window.x().start();
     const int     window_end_x   = window.x().end();
     constexpr int window_step_x  = 16 / sizeof(T);
     Window        window_out     = window;

     // Needed to handle loop left-over
     window_out.set(Window::DimX, Window::Dimension(0, 1, 1));

     switch(pool_info.pool_type)
     {
         case PoolingType::MAX:
             max_poolingMxNxD_fp_neon_ndhwc<T>(src, dst0, pool_info, window_out, window_start_x, window_end_x, window_step_x);
             break;
         case PoolingType::AVG:
             avg_poolingMxNxD_fp_neon_ndhwc<T>(src, dst0, pool_info, window_out, window_start_x, window_end_x, window_step_x);
             break;
         case PoolingType::L2:
             l2_poolingMxNxD_fp_neon_ndhwc<T>(src, dst0, pool_info, window_out, window_start_x, window_end_x, window_step_x);
             break;
         default:
             ARM_COMPUTE_ERROR("Pool operation not supported");
     }
 }

 template void poolingMxNxD_fp_neon_ndhwc<float>(const ITensor *src, ITensor *dst0, Pooling3dLayerInfo &pool_info, const Window &window);
 #if defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(ENABLE_FP16_KERNELS)
 template void poolingMxNxD_fp_neon_ndhwc<float16_t>(const ITensor *src, ITensor *dst0, Pooling3dLayerInfo &pool_info, const Window &window);
 #endif /* defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(ENABLE_FP16_KERNELS) */
 } // namespace cpu
 } // namespace arm_compute
	/*
	* Copyright (c) 2022 Arm Limited.
	*
	* SPDX-License-Identifier: MIT
	*
	* Permission is hereby granted, free of charge, to any person obtaining a copy
	* of this software and associated documentation files (the "Software"), to
	* deal in the Software without restriction, including without limitation the
	* rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
	* sell copies of the Software, and to permit persons to whom the Software is
	* furnished to do so, subject to the following conditions:
	*
	* The above copyright notice and this permission notice shall be included in all
	* copies or substantial portions of the Software.
	*
	* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
	* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
	* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
	* SOFTWARE.
	*/
	#include "arm_compute/core/Helpers.h"
	#include "arm_compute/core/ITensor.h"
	#include "arm_compute/core/Types.h"
	#include "arm_compute/core/utils/misc/Traits.h"
	#include "src/core/NEON/wrapper/intrinsics/intrinsics.h"
	#include "src/core/helpers/WindowHelpers.h"

	#include "src/cpu/kernels/pool3d/neon/impl.h"

	namespace arm_compute
	{
	namespace cpu
	{
	namespace
	{
	inline float calculate_avg_scale(bool exclude_padding, const Coordinates &id, const int pool_size_x, const int pool_size_y, const int pool_size_z, const int upper_bound_w,
	const int upper_bound_h, const int upper_bound_d, const int pad_x, const int pad_y, const int pad_z, const int stride_x, const int stride_y, const int stride_z)
	{
	// Based on NDHWC
	int start_x = id[1] * stride_x - pad_x;
	int start_y = id[2] * stride_y - pad_y;
	int start_z = id[3] * stride_z - pad_z;

	const int end_x = std::min(start_x + pool_size_x, upper_bound_w);
	const int end_y = std::min(start_y + pool_size_y, upper_bound_h);
	const int end_z = std::min(start_z + pool_size_z, upper_bound_d);
	if(exclude_padding)
	{
	start_x = std::max(0, start_x);
	start_y = std::max(0, start_y);
	start_z = std::max(0, start_z);
	}
	return 1.f / ((end_y - start_y) * (end_x - start_x) * (end_z - start_z));
	}


	template <typename T>
	void max_poolingMxNxD_fp_neon_ndhwc(const ITensor src, ITensor dst0, Pooling3dLayerInfo &pool_info, const Window &window_out,
	const int window_start_x, const int window_end_x, const int window_step_x)

	{
	using vtype = wrapper::traits::neon_bitvector<T, wrapper::traits::BitWidth::W128>;
	using vector_type = typename vtype::type;
	using tag_type = typename vtype::tag_type;

	int pool_stride_x = static_cast<int>(pool_info.stride.width);
	int pool_stride_y = static_cast<int>(pool_info.stride.height);
	int pool_stride_z = static_cast<int>(pool_info.stride.depth);

	const int pool_size_x = pool_info.is_global_pooling ? src->info()->tensor_shape().y() : pool_info.pool_size.width;
	const int pool_size_y = pool_info.is_global_pooling ? src->info()->tensor_shape().z() : pool_info.pool_size.height;
	const int pool_size_z = pool_info.is_global_pooling ? src->info()->tensor_shape()[3] : pool_info.pool_size.depth;

	const int pool_pad_top = static_cast<int>(pool_info.padding.top);
	const int pool_pad_left = static_cast<int>(pool_info.padding.left);
	const int pool_pad_front = static_cast<int>(pool_info.padding.front);

	const int input_dim_w = src->info()->dimension(1);
	const int input_dim_h = src->info()->dimension(2);
	const int input_dim_d = src->info()->dimension(3);

	const int y_stride = static_cast<int>(src->info()->strides_in_bytes().y());
	const int z_stride = static_cast<int>(src->info()->strides_in_bytes().z());
	const int w_stride = static_cast<int>(src->info()->strides_in_bytes()[3]);
	const int n_stride = static_cast<int>(src->info()->strides_in_bytes()[4]);

	const uint8_t *in_ptr_start = src->buffer() + src->info()->offset_first_element_in_bytes();

	Iterator out(dst0, window_out);

	vector_type vres;
	execute_window_loop(window_out, [&](const Coordinates & id)
	{
	// Computing the theoretical input starting/ending points
	const int in_idx_width = static_cast<int>(id.y()) * pool_stride_x - pool_pad_left;
	const int in_idx_height = static_cast<int>(id.z()) * pool_stride_y - pool_pad_top;
	const int in_idx_depth = static_cast<int>(id[3]) * pool_stride_z - pool_pad_front;

	const int pool_start_x = std::max(0, -in_idx_width);
	const int pool_end_x_t = std::min(input_dim_w + pool_pad_left - in_idx_width, pool_size_x);
	const int pool_start_y = std::max(0, -in_idx_height);
	const int pool_end_y_t = std::min(input_dim_h + pool_pad_top - in_idx_height, pool_size_y);

	const int pool_start_z = std::max(0, -in_idx_depth);
	const int pool_end_z_t = std::min(input_dim_d + pool_pad_front - in_idx_depth, pool_size_z);

	// The end of width to consider in calculation should exclude PAD_X, PAD_Y and PAD_Z
	const int pool_end_x = std::min(pool_end_x_t, input_dim_w - in_idx_width);
	const int pool_end_y = std::min(pool_end_y_t, input_dim_h - in_idx_height);
	const int pool_end_z = std::min(pool_end_z_t, input_dim_d - in_idx_depth);

	const uint8_t in_ptr_n = in_ptr_start + id[4] n_stride;

	int x_off = window_start_x;

	for(; x_off <= (window_end_x - window_step_x); x_off += window_step_x) // C
	{
	vres = wrapper::vdup_n(static_cast<T>(-std::numeric_limits<float>::infinity()), tag_type());
	for(int z = pool_start_z; z < pool_end_z; ++z)
	{
	const uint8_t in_ptr_z = in_ptr_n + (z + in_idx_depth) w_stride;
	for(int y = pool_start_y; y < pool_end_y; ++y)
	{
	const uint8_t in_ptr_y = in_ptr_z + (y + in_idx_height) z_stride;
	for(int x = pool_start_x; x < pool_end_x; ++x)
	{
	const uint8_t in_ptr_x = in_ptr_y + (x + in_idx_width) y_stride;
	const vector_type data = wrapper::vloadq(reinterpret_cast<const T *>(in_ptr_x) + x_off);
	vres = wrapper::vmax(vres, data);
	}
	}
	}
	// Store result
	wrapper::vstore(reinterpret_cast<T *>(out.ptr()) + x_off, vres);
	}

	// Left-overs loop
	for(; x_off < window_end_x; ++x_off)
	{
	T res(0);
	res = -std::numeric_limits<float>::infinity();
	for(int z = pool_start_z; z < pool_end_z; ++z)
	{
	const uint8_t in_ptr_z = in_ptr_n + (z + in_idx_depth) w_stride;
	for(int y = pool_start_y; y < pool_end_y; ++y)
	{
	const uint8_t in_ptr_y = in_ptr_z + (y + in_idx_height) z_stride;
	for(int x = pool_start_x; x < pool_end_x; ++x)
	{
	const uint8_t in_ptr_x = in_ptr_y + (x + in_idx_width) y_stride;
	const T data = (reinterpret_cast<const T >(in_ptr_x) + x_off);
	res = std::max(res, data);
	}
	}
	}
	// Store result
	(reinterpret_cast<T >(out.ptr()) + x_off) = res;
	}
	},
	out);
	}

	template <typename T>
	void avg_poolingMxNxD_fp_neon_ndhwc(const ITensor src, ITensor dst0, Pooling3dLayerInfo &pool_info,
	const Window &window_out, const int window_start_x, const int window_end_x, const int window_step_x)
	{
	using vtype = wrapper::traits::neon_bitvector<T, wrapper::traits::BitWidth::W128>;
	using vector_type = typename vtype::type;
	using tag_type = typename vtype::tag_type;

	int pool_stride_x = static_cast<int>(pool_info.stride.width);
	int pool_stride_y = static_cast<int>(pool_info.stride.height);
	int pool_stride_z = static_cast<int>(pool_info.stride.depth);

	const int pool_size_x = pool_info.is_global_pooling ? src->info()->tensor_shape().y() : pool_info.pool_size.width;
	const int pool_size_y = pool_info.is_global_pooling ? src->info()->tensor_shape().z() : pool_info.pool_size.height;
	const int pool_size_z = pool_info.is_global_pooling ? src->info()->tensor_shape()[3] : pool_info.pool_size.depth;

	const int pool_pad_top = static_cast<int>(pool_info.padding.top);
	const int pool_pad_bottom = static_cast<int>(pool_info.padding.bottom);
	const int pool_pad_left = static_cast<int>(pool_info.padding.left);
	const int pool_pad_right = static_cast<int>(pool_info.padding.right);
	const int pool_pad_front = static_cast<int>(pool_info.padding.front);
	const int pool_pad_back = static_cast<int>(pool_info.padding.back);

	const int upper_bound_w = src->info()->dimension(1) + (pool_info.exclude_padding ? 0 : pool_pad_right);
	const int upper_bound_h = src->info()->dimension(2) + (pool_info.exclude_padding ? 0 : pool_pad_bottom);
	const int upper_bound_d = src->info()->dimension(3) + (pool_info.exclude_padding ? 0 : pool_pad_back);

	const int input_dim_w = src->info()->dimension(1);
	const int input_dim_h = src->info()->dimension(2);
	const int input_dim_d = src->info()->dimension(3);

	const int y_stride = static_cast<int>(src->info()->strides_in_bytes().y());
	const int z_stride = static_cast<int>(src->info()->strides_in_bytes().z());
	const int w_stride = static_cast<int>(src->info()->strides_in_bytes()[3]);
	const int n_stride = static_cast<int>(src->info()->strides_in_bytes()[4]);

	const uint8_t *in_ptr_start = src->buffer() + src->info()->offset_first_element_in_bytes();

	Iterator out(dst0, window_out);

	vector_type vres;
	execute_window_loop(window_out, [&](const Coordinates & id)
	{
	// Computing the theoretical input starting/ending points
	const int in_idx_width = static_cast<int>(id.y()) * pool_stride_x - pool_pad_left;
	const int in_idx_height = static_cast<int>(id.z()) * pool_stride_y - pool_pad_top;
	const int in_idx_depth = static_cast<int>(id[3]) * pool_stride_z - pool_pad_front;

	const int pool_start_x = std::max(0, -in_idx_width);
	const int pool_end_x_t = std::min(input_dim_w + pool_pad_left - in_idx_width, pool_size_x);
	const int pool_start_y = std::max(0, -in_idx_height);
	const int pool_end_y_t = std::min(input_dim_h + pool_pad_top - in_idx_height, pool_size_y);

	const int pool_start_z = std::max(0, -in_idx_depth);
	const int pool_end_z_t = std::min(input_dim_d + pool_pad_front - in_idx_depth, pool_size_z);

	// The end of width to consider in calculation should exclude PAD_X, PAD_Y and PAD_Z
	const int pool_end_x = std::min(pool_end_x_t, input_dim_w - in_idx_width);
	const int pool_end_y = std::min(pool_end_y_t, input_dim_h - in_idx_height);
	const int pool_end_z = std::min(pool_end_z_t, input_dim_d - in_idx_depth);

	const uint8_t in_ptr_n = in_ptr_start + id[4] n_stride;

	// Calculate scale
	const float scale = calculate_avg_scale(pool_info.exclude_padding, id, pool_size_x, pool_size_y, pool_size_z, upper_bound_w, upper_bound_h, upper_bound_d, pool_pad_left,
	pool_pad_top, pool_pad_front, pool_stride_x,
	pool_stride_y, pool_stride_z);
	const vector_type scale_v = wrapper::vdup_n(static_cast<T>(scale), tag_type());

	int x_off = window_start_x;

	for(; x_off <= (window_end_x - window_step_x); x_off += window_step_x) // C
	{
	// Perform pooling
	vres = wrapper::vdup_n(static_cast<T>(0.0f), tag_type());
	for(int z = pool_start_z; z < pool_end_z; ++z)
	{
	const uint8_t in_ptr_z = in_ptr_n + (z + in_idx_depth) w_stride;
	for(int y = pool_start_y; y < pool_end_y; ++y)
	{
	const uint8_t in_ptr_y = in_ptr_z + (y + in_idx_height) z_stride;
	for(int x = pool_start_x; x < pool_end_x; ++x)
	{
	const uint8_t in_ptr_x = in_ptr_y + (x + in_idx_width) y_stride;
	const vector_type data = wrapper::vloadq(reinterpret_cast<const T *>(in_ptr_x) + x_off);
	vres = wrapper::vadd(vres, data);
	}
	}
	}

	// Divide by scale
	vres = wrapper::vmul(vres, scale_v);

	// Store result
	wrapper::vstore(reinterpret_cast<T *>(out.ptr()) + x_off, vres);
	}

	// Left-overs loop
	for(; x_off < window_end_x; ++x_off)
	{
	T res(0);

	for(int z = pool_start_z; z < pool_end_z; ++z)
	{
	const uint8_t in_ptr_z = in_ptr_n + (z + in_idx_depth) w_stride;
	for(int y = pool_start_y; y < pool_end_y; ++y)
	{
	const uint8_t in_ptr_y = in_ptr_z + (y + in_idx_height) z_stride;
	for(int x = pool_start_x; x < pool_end_x; ++x)
	{
	const uint8_t in_ptr_x = in_ptr_y + (x + in_idx_width) y_stride;
	const T data = (reinterpret_cast<const T >(in_ptr_x) + x_off);
	res += data;
	}
	}
	}

	// Divide by scale
	res *= scale;

	// Store result
	(reinterpret_cast<T >(out.ptr()) + x_off) = res;
	}
	},
	out);
	}

	template <typename T>
	void l2_poolingMxNxD_fp_neon_ndhwc(const ITensor src, ITensor dst0, Pooling3dLayerInfo &pool_info,
	const Window &window_out, const int window_start_x, const int window_end_x, const int window_step_x)
	{
	using vtype = wrapper::traits::neon_bitvector<T, wrapper::traits::BitWidth::W128>;
	using vector_type = typename vtype::type;
	using tag_type = typename vtype::tag_type;

	int pool_stride_x = static_cast<int>(pool_info.stride.width);
	int pool_stride_y = static_cast<int>(pool_info.stride.height);
	int pool_stride_z = static_cast<int>(pool_info.stride.depth);

	const int pool_size_x = pool_info.is_global_pooling ? src->info()->tensor_shape().y() : pool_info.pool_size.width;
	const int pool_size_y = pool_info.is_global_pooling ? src->info()->tensor_shape().z() : pool_info.pool_size.height;
	const int pool_size_z = pool_info.is_global_pooling ? src->info()->tensor_shape()[3] : pool_info.pool_size.depth;

	const int pool_pad_top = static_cast<int>(pool_info.padding.top);
	const int pool_pad_bottom = static_cast<int>(pool_info.padding.bottom);
	const int pool_pad_left = static_cast<int>(pool_info.padding.left);
	const int pool_pad_right = static_cast<int>(pool_info.padding.right);
	const int pool_pad_front = static_cast<int>(pool_info.padding.front);
	const int pool_pad_back = static_cast<int>(pool_info.padding.back);

	const int upper_bound_w = src->info()->dimension(1) + (pool_info.exclude_padding ? 0 : pool_pad_right);
	const int upper_bound_h = src->info()->dimension(2) + (pool_info.exclude_padding ? 0 : pool_pad_bottom);
	const int upper_bound_d = src->info()->dimension(3) + (pool_info.exclude_padding ? 0 : pool_pad_back);

	const int input_dim_w = src->info()->dimension(1);
	const int input_dim_h = src->info()->dimension(2);
	const int input_dim_d = src->info()->dimension(3);

	const int y_stride = static_cast<int>(src->info()->strides_in_bytes().y());
	const int z_stride = static_cast<int>(src->info()->strides_in_bytes().z());
	const int w_stride = static_cast<int>(src->info()->strides_in_bytes()[3]);
	const int n_stride = static_cast<int>(src->info()->strides_in_bytes()[4]);

	const uint8_t *in_ptr_start = src->buffer() + src->info()->offset_first_element_in_bytes();

	Iterator out(dst0, window_out);

	vector_type vres;
	execute_window_loop(window_out, [&](const Coordinates & id)
	{
	// Computing the theoretical input starting/ending points
	const int in_idx_width = static_cast<int>(id.y()) * pool_stride_x - pool_pad_left;
	const int in_idx_height = static_cast<int>(id.z()) * pool_stride_y - pool_pad_top;
	const int in_idx_depth = static_cast<int>(id[3]) * pool_stride_z - pool_pad_front;

	const int pool_start_x = std::max(0, -in_idx_width);
	const int pool_end_x_t = std::min(input_dim_w + pool_pad_left - in_idx_width, pool_size_x);
	const int pool_start_y = std::max(0, -in_idx_height);
	const int pool_end_y_t = std::min(input_dim_h + pool_pad_top - in_idx_height, pool_size_y);

	const int pool_start_z = std::max(0, -in_idx_depth);
	const int pool_end_z_t = std::min(input_dim_d + pool_pad_front - in_idx_depth, pool_size_z);

	// The end of width to consider in calculation should exclude PAD_X, PAD_Y and PAD_Z
	const int pool_end_x = std::min(pool_end_x_t, input_dim_w - in_idx_width);
	const int pool_end_y = std::min(pool_end_y_t, input_dim_h - in_idx_height);
	const int pool_end_z = std::min(pool_end_z_t, input_dim_d - in_idx_depth);

	const uint8_t in_ptr_n = in_ptr_start + id[4] n_stride;

	// Calculate scale
	const float scale = calculate_avg_scale(pool_info.exclude_padding, id, pool_size_x, pool_size_y, pool_size_z, upper_bound_w, upper_bound_h, upper_bound_d, pool_pad_left,
	pool_pad_top, pool_pad_front, pool_stride_x,
	pool_stride_y, pool_stride_z);

	int x_off = window_start_x;

	for(; x_off <= (window_end_x - window_step_x); x_off += window_step_x) // C
	{
	// Perform pooling
	vres = wrapper::vdup_n(static_cast<T>(0.0f), tag_type());
	for(int z = pool_start_z; z < pool_end_z; ++z)
	{
	const uint8_t in_ptr_z = in_ptr_n + (z + in_idx_depth) w_stride;
	for(int y = pool_start_y; y < pool_end_y; ++y)
	{
	const uint8_t in_ptr_y = in_ptr_z + (y + in_idx_height) z_stride;
	for(int x = pool_start_x; x < pool_end_x; ++x)
	{
	const uint8_t in_ptr_x = in_ptr_y + (x + in_idx_width) y_stride;
	const vector_type data = wrapper::vloadq(reinterpret_cast<const T *>(in_ptr_x) + x_off);
	vres = wrapper::vmla(vres, data, data);
	}
	}
	}

	const vector_type scale_v = wrapper::vdup_n(static_cast<T>(scale), tag_type());

	// Divide by scale
	vres = wrapper::vmul(vres, scale_v);

	// Calculate square-root
	vres = wrapper::vinv(wrapper::vinvsqrt(vres));

	// Store result
	wrapper::vstore(reinterpret_cast<T *>(out.ptr()) + x_off, vres);
	}

	// Left-overs loop
	for(; x_off < window_end_x; ++x_off)
	{
	T res(0);

	for(int z = pool_start_z; z < pool_end_z; ++z)
	{
	const uint8_t in_ptr_z = in_ptr_n + (z + in_idx_depth) w_stride;
	for(int y = pool_start_y; y < pool_end_y; ++y)
	{
	const uint8_t in_ptr_y = in_ptr_z + (y + in_idx_height) z_stride;
	for(int x = pool_start_x; x < pool_end_x; ++x)
	{
	const uint8_t in_ptr_x = in_ptr_y + (x + in_idx_width) y_stride;
	const T data = (reinterpret_cast<const T >(in_ptr_x) + x_off);
	res += data * data;
	}
	}
	}

	// Divide by scale
	res *= scale;

	// Square root
	res = std::sqrt(res);

	// Store result
	(reinterpret_cast<T >(out.ptr()) + x_off) = res;
	}
	},
	out);
	}
	} // namespace

	template <typename T>
	void poolingMxNxD_fp_neon_ndhwc(const ITensor src, ITensor dst0, Pooling3dLayerInfo &pool_info, const Window &window)
	{
	const int window_start_x = window.x().start();
	const int window_end_x = window.x().end();
	constexpr int window_step_x = 16 / sizeof(T);
	Window window_out = window;

	// Needed to handle loop left-over
	window_out.set(Window::DimX, Window::Dimension(0, 1, 1));

	switch(pool_info.pool_type)
	{
	case PoolingType::MAX:
	max_poolingMxNxD_fp_neon_ndhwc<T>(src, dst0, pool_info, window_out, window_start_x, window_end_x, window_step_x);
	break;
	case PoolingType::AVG:
	avg_poolingMxNxD_fp_neon_ndhwc<T>(src, dst0, pool_info, window_out, window_start_x, window_end_x, window_step_x);
	break;
	case PoolingType::L2:
	l2_poolingMxNxD_fp_neon_ndhwc<T>(src, dst0, pool_info, window_out, window_start_x, window_end_x, window_step_x);
	break;
	default:
	ARM_COMPUTE_ERROR("Pool operation not supported");
	}
	}

	template void poolingMxNxD_fp_neon_ndhwc<float>(const ITensor src, ITensor dst0, Pooling3dLayerInfo &pool_info, const Window &window);
	#if defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(ENABLE_FP16_KERNELS)
	template void poolingMxNxD_fp_neon_ndhwc<float16_t>(const ITensor src, ITensor dst0, Pooling3dLayerInfo &pool_info, const Window &window);
	#endif /* defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(ENABLE_FP16_KERNELS) */
	} // namespace cpu
	} // namespace arm_compute