text: #pragma once
  
  #if THRUST_DEVICE_COMPILER == THRUST_DEVICE_COMPILER_NVCC
  #include <thrust/system/cuda/config.h>
  
  #include <thrust/detail/cstdint.h>
  #include <thrust/detail/temporary_array.h>
  #include <thrust/detail/raw_reference_cast.h>
  #include <thrust/detail/type_traits/iterator/is_output_iterator.h>
  #include <cub/device/device_reduce.cuh>
  #include <thrust/functional.h>
  #include <thrust/system/ktt/ktt_functors.h>
  #include <thrust/detail/minmax.h>
  #include <thrust/distance.h>
  #include <thrust/detail/alignment.h>
  #include <thrust/device_vector.h>
  
  #include <Ktt.h>
  
  namespace thrust
  {
  namespace cuda_ktt
  {
  
  const std::string reductionSource = std::string("\n") +
  "// todo: float to T\n" +
  "/*#if VECTOR_SIZE == 1\n" +
  "typedef float VEC;\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 2\n" +
  "typedef float2 VEC;\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 4\n" +
  "typedef float4 VEC;\n" +
  "#endif*/\n" +
  "\n" +
  "// todo: float to T\n" +
  "inline __host__ __device__ void operator+=(float2 &a, float2 b)\n" +
  "{\n" +
  "    a.x += b.x;\n" +
  "    a.y += b.y;\n" +
  "}\n" +
  "\n" +
  "inline __host__ __device__ void operator+=(float4 &a, float4 b)\n" +
  "{\n" +
  "    a.x += b.x;\n" +
  "    a.y += b.y;\n" +
  "    a.z += b.z;\n" +
  "    a.w += b.w;\n" +
  "}\n" +
  "\n" +
  "template <typename T, typename Operator>\n" +
  "__global__ void reduce(const T* in, T* out, unsigned int n, unsigned int inOffset, unsigned int outOffset)\n" +
  "{\n" +
  "    Operator op;\n" +
  "    unsigned int tid = threadIdx.x;\n" +
  "    unsigned int i = blockDim.x*blockIdx.x + threadIdx.x;\n" +
  "\n" +
  "    __shared__ float buf[WORK_GROUP_SIZE_X];\n" +
  "\n" +
  "#if UNBOUNDED_WG == 1\n" +
  "    unsigned int addr = i + inOffset;\n" +
  "\n" +
  "    if (i < (n + VECTOR_SIZE - 1) / VECTOR_SIZE)\n" +
  "    {\n" +
  "#if VECTOR_SIZE == 1\n" +
  "        buf[tid] = in[addr];\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 2\n" +
  "        buf[tid] = op(in[addr].x, in[addr].y);\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 4\n" +
  "        buf[tid] = op(in[addr].x, op(in[addr].y, op(in[addr].z, in[addr].w)));\n" +
  "#endif\n" +
  "    }\n" +
  "    else\n" +
  "    {\n" +
  "        buf[tid] = 0.0f;\n" +
  "    }\n" +
  "#else /*UNBOUNDED_WG != 1*/\n" +
  "    T partial;\n" +
  "#if VECTOR_SIZE == 1\n" +
  "    partial = 0.0f;\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 2\n" +
  "    partial.x = 0.0f;\n" +
  "    partial.y = 0.0f;\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 4\n" +
  "    partial.x = 0.0f;\n" +
  "    partial.y = 0.0f;\n" +
  "    partial.z = 0.0f;\n" +
  "    partial.w = 0.0f;\n" +
  "#endif\n" +
  "\n" +
  "    while (i < n / VECTOR_SIZE)\n" +
  "    {\n" +
  "        partial = op(partial, in[i + inOffset]);\n" +
  "        i += WORK_GROUP_SIZE_X * WG_NUM;\n" +
  "    }\n" +
  "#if VECTOR_SIZE == 1\n" +
  "    buf[tid] = partial;\n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 2\n" +
  "    buf[tid] = op(partial.x, partial.y);    \n" +
  "#endif\n" +
  "#if VECTOR_SIZE == 4\n" +
  "    buf[tid] = op(partial.x, op(partial.y, op(partial.z, partial.w)));\n" +
  "#endif\n" +
  "#endif /* UNBOUNDED_WG != 1 */\n" +
  "\n" +
  "    __syncthreads();\n" +
  "\n" +
  "#if WORK_GROUP_SIZE_X >= 512\n" +
  "    if (tid < 256)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 256]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 256\n" +
  "    if (tid < 128)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 128]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 128\n" +
  "    if (tid < 64)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 64]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 64\n" +
  "    if (tid < 32)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 32]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 32\n" +
  "    if (tid < 16)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 16]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 16\n" +
  "    if (tid < 8)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 8]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 8\n" +
  "    if (tid < 4)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 4]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 4\n" +
  "    if (tid < 2)\n" +
  "        buf[tid] = op(buf[tid], buf[tid + 2]);\n" +
  "    __syncthreads();\n" +
  "#endif\n" +
  "#if WORK_GROUP_SIZE_X >= 2\n" +
  "    if (tid < 1)\n" +
  "    {\n" +
  "        buf[0] = op(buf[0], buf[1]);\n" +
  "    }\n" +
  "#endif\n" +
  "#if USE_ATOMICS == 1\n" +
  "    if (tid < 1)\n" +
  "        // todo: figure out how to make this compatible with generic operator\n" +
  "        atomicAdd(out + outOffset, buf[0]);\n" +
  "#else\n" +
  "    // the last group\n" +
  "    if (blockIdx.x == gridDim.x-1)\n" +
  "    {\n" +
  "        if (tid == 0)\n" +
  "            out[blockIdx.x + outOffset] = buf[0];\n" +
  "        else if (tid < VECTOR_SIZE)\n" +
  "            out[blockIdx.x + outOffset + tid] = 0.0f;\n" +
  "    }\n" +
  "    // other groups\n" +
  "    else\n" +
  "    {\n" +
  "        if (tid == 0)\n" +
  "            out[blockIdx.x + outOffset] = buf[0];\n" +
  "    }\n" +
  "#endif\n" +
  "}";
  
  template <typename T>
  struct ReductionData
  {
      ReductionData() : m_IsInitialized(false)
      {}
  
      void SetArguments()
      {
          auto& tuner = GetTuner();
          tuner.SetArguments(m_ReductionDefinition, {m_SrcId, m_DstId, m_NId, m_InOffsetId, m_OutOffsetId});
      }
  
      void ClearArguments()
      {
          auto& tuner = GetTuner();
          tuner.SetArguments(m_ReductionDefinition, {});
      }
  
      ktt::KernelDefinitionId m_ReductionDefinition;
      ktt::KernelId m_ReductionKernel;
      ktt::ArgumentId m_SrcId;
      ktt::ArgumentId m_DstId;
      ktt::ArgumentId m_NId;
      ktt::ArgumentId m_InOffsetId;
      ktt::ArgumentId m_OutOffsetId;
      device_vector<T> m_Dst;
      std::string m_OpSource;
      bool m_IsInitialized;
  };
  
  template <typename T>
  ReductionData<T>& GetReductionData()
  {
      static ReductionData<T> data;
      return data;
  }
  
  ktt::Tuner& GetTuner()
  {
      static ktt::Tuner tuner(0, 0, ktt::ComputeApi::CUDA);
      return tuner;
  }
  
  template <typename InputIt, typename Size, typename T, typename BinaryOp>
  void InitializeReduction(InputIt first, Size numItems, BinaryOp binaryOp)
  {
      auto& reductionData = GetReductionData<T>();
      auto& tuner = GetTuner();
  
      if (reductionData.m_IsInitialized)
      {
          reductionData.ClearArguments();
          tuner.RemoveArgument(reductionData.m_SrcId);
          T* srcBuffer = thrust::raw_pointer_cast(&(*first));
          const size_t size = numItems * sizeof(T);
          reductionData.m_SrcId = tuner.AddArgumentVector<T>(static_cast<ktt::ComputeBuffer>(srcBuffer), size, ktt::ArgumentAccessType::ReadWrite, ktt::ArgumentMemoryLocation::Device);
          reductionData.SetArguments();
  
          if (reductionData.m_Dst.size() < numItems)
          {
              reductionData.m_Dst.resize(numItems, static_cast<T>(0));
          }
  
          return;
      }
  
      const size_t size = numItems * sizeof(T);
      const uint32_t nUp = ((numItems + 512 - 1) / 512) * 512; // maximum WG size used in tuning parameters
      const uint32_t nAlloc = ((numItems + 16 - 1) / 16) * 16; // pad to the longest vector size
  
      tuner.SetGlobalSizeType(ktt::GlobalSizeType::OpenCL);
      tuner.SetTimeUnit(ktt::TimeUnit::Microseconds);
  
      //static_assert(IsCudaFunctor<BinaryOp>::m_Value, "test failed");
      reductionData.m_OpSource = binaryOp.GetCudaCode();
      reductionData.m_ReductionDefinition = tuner.AddKernelDefinition("reduce", reductionData.m_OpSource + reductionSource, ktt::DimensionVector(nUp), ktt::DimensionVector(), {"float", "plus_ktt<float>"});
      reductionData.m_ReductionKernel = tuner.CreateSimpleKernel("Reduction", reductionData.m_ReductionDefinition);
  
      T* srcBuffer = thrust::raw_pointer_cast(&(*first));
      reductionData.m_SrcId = tuner.AddArgumentVector<T>(static_cast<ktt::ComputeBuffer>(srcBuffer), size, ktt::ArgumentAccessType::ReadWrite, ktt::ArgumentMemoryLocation::Device);
  
      auto& dst = reductionData.m_Dst;
      dst.resize(numItems, static_cast<T>(0));
      T* dstBuffer = thrust::raw_pointer_cast(dst.data());
      reductionData.m_DstId = tuner.AddArgumentVector<T>(static_cast<ktt::ComputeBuffer>(dstBuffer), size, ktt::ArgumentAccessType::ReadWrite, ktt::ArgumentMemoryLocation::Device);
  
      reductionData.m_NId = tuner.AddArgumentScalar(numItems);
      uint32_t offset = 0;
      reductionData.m_InOffsetId = tuner.AddArgumentScalar(offset);
      reductionData.m_OutOffsetId = tuner.AddArgumentScalar(offset);
      reductionData.SetArguments();
  
      // get number of compute units
      const ktt::DeviceInfo di = tuner.GetCurrentDeviceInfo();
      //std::cout << "Number of compute units: " << di.GetMaxComputeUnits() << std::endl;
      size_t cus = di.GetMaxComputeUnits();
  
      tuner.AddParameter(reductionData.m_ReductionKernel, "WORK_GROUP_SIZE_X", std::vector<uint64_t>{32, 64, 96, 128, 160, 196, 224, 256, 288, 320, 352, 384, 416, 448, 480, 512, 576, 640, 704, 768, 832, 896, 960, 1024});
      tuner.AddThreadModifier(reductionData.m_ReductionKernel, {reductionData.m_ReductionDefinition}, ktt::ModifierType::Local, ktt::ModifierDimension::X, "WORK_GROUP_SIZE_X", ktt::ModifierAction::Multiply);
      tuner.AddParameter(reductionData.m_ReductionKernel, "UNBOUNDED_WG", std::vector<uint64_t>{0, 1});
      tuner.AddParameter(reductionData.m_ReductionKernel, "WG_NUM", std::vector<uint64_t>{0, cus, cus * 2, cus * 3, cus * 4, cus * 5, cus * 6, cus * 7, cus * 8, cus * 10, cus * 12, cus * 14, cus * 16, cus * 20, cus * 24, cus * 28, cus * 32, cus * 40, cus * 48, cus * 56, cus * 64});
  
      tuner.AddParameter(reductionData.m_ReductionKernel, "VECTOR_SIZE", std::vector<uint64_t>{1/*, 2, 4*/});
  
      tuner.AddThreadModifier(reductionData.m_ReductionKernel, {reductionData.m_ReductionDefinition}, ktt::ModifierType::Global, ktt::ModifierDimension::X, "VECTOR_SIZE", ktt::ModifierAction::Divide);
      tuner.AddParameter(reductionData.m_ReductionKernel, "USE_ATOMICS", std::vector<uint64_t>{0/*, 1*/});
  
      auto persistConstraint = [](const std::vector<uint64_t>& v) {return (v[0] && v[1] == 0) || (!v[0] && v[1] > 0); };
      tuner.AddConstraint(reductionData.m_ReductionKernel, {"UNBOUNDED_WG", "WG_NUM"}, persistConstraint);
      auto persistentAtomic = [](const std::vector<uint64_t>& v) {return (v[0] == 1) || (v[0] == 0 && v[1] == 1); };
      tuner.AddConstraint(reductionData.m_ReductionKernel, {"UNBOUNDED_WG", "USE_ATOMICS"}, persistentAtomic);
      auto unboundedWG = [](const std::vector<uint64_t>& v) {return (!v[0] || v[1] >= 32); };
      tuner.AddConstraint(reductionData.m_ReductionKernel, {"UNBOUNDED_WG", "WORK_GROUP_SIZE_X"}, unboundedWG);
  
      /*if constexpr (!rapidTest)
      {
      tuner.SetReferenceComputation(dstId, [&src](void* buffer)
      {
      float* result = static_cast<float*>(buffer);
      std::vector<double> resD(src.size());
      size_t resSize = src.size();
  
      for (size_t i = 0; i < resSize; ++i)
      {
      resD[i] = static_cast<double>(src[i]);
      }
  
      while (resSize > 1)
      {
      for (size_t i = 0; i < resSize / 2; ++i)
      {
      resD[i] = resD[i * 2] + resD[i * 2 + 1];
      }
  
      if (resSize % 2 != 0)
      {
      resD[resSize / 2 - 1] += resD[resSize - 1];
      }
  
      resSize = resSize / 2;
      }
  
      std::cout << "Reference in double: " << std::setprecision(10) << resD[0] << std::endl;
      result[0] = static_cast<float>(resD[0]);
      });
  
      tuner.SetValidationMethod(ktt::ValidationMethod::SideBySideComparison, static_cast<double>(n) * 10'000.0 / 10'000'000.0);
      tuner.SetValidationRange(dstId, 1);
      }*/
  
      tuner.SetLauncher(reductionData.m_ReductionKernel, [&reductionData](ktt::ComputeInterface& interface)
      {
          const ktt::DimensionVector& globalSize = interface.GetCurrentGlobalSize(reductionData.m_ReductionDefinition);
          const ktt::DimensionVector& localSize = interface.GetCurrentLocalSize(reductionData.m_ReductionDefinition);
          const std::vector<ktt::ParameterPair>& pairs = interface.GetCurrentConfiguration().GetPairs();
          ktt::DimensionVector myGlobalSize = globalSize;
  
          // change global size for constant numbers of work-groups
          // this may be done by thread modifier operators as well
          if (ktt::ParameterPair::GetParameterValue<uint64_t>(pairs, "UNBOUNDED_WG") == 0)
          {
              myGlobalSize = ktt::DimensionVector(ktt::ParameterPair::GetParameterValue<uint64_t>(pairs, "WG_NUM")
                  * localSize.GetSizeX());
          }
  
          interface.RunKernel(reductionData.m_ReductionDefinition, myGlobalSize, localSize);
  
          // execute kernel log n times, when atomics are not used 
          if (ktt::ParameterPair::GetParameterValue<uint64_t>(pairs, "USE_ATOMICS") == 0)
          {
              uint32_t n = static_cast<uint32_t>(globalSize.GetSizeX() / localSize.GetSizeX());
              uint32_t inOffset = 0;
              uint32_t outOffset = n;
              uint32_t vectorSize = static_cast<uint32_t>(ktt::ParameterPair::GetParameterValue<uint64_t>(pairs, "VECTOR_SIZE"));
              uint32_t wgSize = static_cast<uint32_t>(localSize.GetSizeX());
              size_t iterations = 0; // make sure the end result is in the correct buffer
  
              while (n > 1 || iterations % 2 == 1)
              {
                  interface.SwapArguments(reductionData.m_ReductionDefinition, reductionData.m_SrcId, reductionData.m_DstId);
                  myGlobalSize.SetSizeX((n + vectorSize - 1) / vectorSize);
                  myGlobalSize.SetSizeX(((myGlobalSize.GetSizeX() - 1) / wgSize + 1) * wgSize);
  
                  if (myGlobalSize == localSize)
                  {
                      outOffset = 0; // only one WG will be executed
                  }
  
                  interface.UpdateScalarArgument(reductionData.m_NId, &n);
                  interface.UpdateScalarArgument(reductionData.m_OutOffsetId, &outOffset);
                  interface.UpdateScalarArgument(reductionData.m_InOffsetId, &inOffset);
  
                  interface.RunKernel(reductionData.m_ReductionDefinition, myGlobalSize, localSize);
                  n = (n + wgSize * vectorSize - 1) / (wgSize * vectorSize);
                  inOffset = outOffset / vectorSize; // input is vectorized, output is scalar
                  outOffset += n;
                  ++iterations;
              }
          }
      });
  
      reductionData.m_IsInitialized = true;
  }
  
  __thrust_exec_check_disable__
  template <typename Derived, typename InputIt, typename Size, typename T, typename BinaryOp>
  __host__ __device__ T reduce_n(execution_policy<Derived>& policy, InputIt first, Size num_items, T init, BinaryOp binary_op)
  {
      if (__THRUST_HAS_CUDART__)
          return reduce_n_impl(policy, first, num_items, init, binary_op);
  
  #if !__THRUST_HAS_CUDART__
      return reduce(cvt_to_seq(derived_cast(policy)), first, first + num_items, init, binary_op);
  #endif
  }
  
  template <typename Derived,
            typename InputIt,
            typename Size,
            typename T,
            typename BinaryOp>
  T reduce_n_impl(execution_policy<Derived>& policy,
                  InputIt                    first,
                  Size                       num_items,
                  T                          init,
                  BinaryOp                   binary_op)
  {
      InitializeReduction<InputIt, Size, T, BinaryOp>(first, num_items, binary_op);
      auto& tuner = GetTuner();
      auto& reductionData = GetReductionData<T>();
  
      T result = 0;
      ktt::BufferOutputDescriptor output(reductionData.m_DstId, &result, sizeof(T));
      const auto results = tuner.TuneIteration(reductionData.m_ReductionKernel, {output});
      return result;
  }
  
  } // end namespace cuda_ktt
  } // end namespace thrust
  
  #include <thrust/memory.h>
  #include <thrust/reduce.h>
  
  #endif
  