                                                       THCudaTensor_size(state, inputImages, 3),
                                                       THCudaTensor_size(state, inputImages, 1),
                                                       THCudaTensor_size(state, inputImages, 2),
+                                                      THCudaTensor_size(state, output, 1),
                                                       THCudaTensor_size(state, output, 2));
 
 

   // threads : used for features

   const int xOut = blockIdx.x*blockDim.y+threadIdx.y;
   const bool withinImageBounds = xOut < gradOutput_width;
   const bool withinGridBounds = blockIdx.x*blockDim.y + threadIdx.x / 2 < gradOutput_width;

   const int yOut = blockIdx.y;
   const int width = inputImages_width;
   const int height = inputImages_height;

   const int b = blockIdx.z;

   float yf,xf;

   __shared__ float gridData[32];
   if (threadIdx.y==0 && withinGridBounds)
   {
      gridData[threadIdx.x] = grids_data[b*grids_strideBatch + yOut*grids_strideHeight + xOut*grids_strideWidth + threadIdx.x];
   }
   __syncthreads();

   if(withinImageBounds)
   {
      yf = gridData[threadIdx.y*2];
      xf = gridData[threadIdx.y*2+1];

      int yInTopLeft, xInTopLeft;
      float yWeightTopLeft, xWeightTopLeft;
      getTopLeft(xf, inputImages_width, xInTopLeft, xWeightTopLeft);
      getTopLeft(yf, inputImages_height, yInTopLeft, yWeightTopLeft);

      xWeightTopLeft = 0.5;
      yWeightTopLeft = 0.5;

      const int inTopLeftAddress = inputImages_strideBatch * b + inputImages_strideHeight * yInTopLeft + inputImages_strideWidth * xInTopLeft;
      const int inTopRightAddress = inTopLeftAddress + inputImages_strideWidth;
      const int inBottomLeftAddress = inTopLeftAddress + inputImages_strideHeight;
      const int inBottomRightAddress = inBottomLeftAddress + inputImages_strideWidth;

      const int gradInputImagesTopLeftAddress = gradInputImages_strideBatch * b + gradInputImages_strideHeight * yInTopLeft + gradInputImages_strideWidth * xInTopLeft;
      const int gradInputImagesTopRightAddress = gradInputImagesTopLeftAddress + gradInputImages_strideWidth;
      const int gradInputImagesBottomLeftAddress = gradInputImagesTopLeftAddress + gradInputImages_strideHeight;
      const int gradInputImagesBottomRightAddress = gradInputImagesBottomLeftAddress + gradInputImages_strideWidth;

      const int inTopLeftMaskAddress = masks_strideBatch * b + masks_strideHeight * yInTopLeft + masks_strideWidth * xInTopLeft;
      const int inTopRightMaskAddress = inTopLeftMaskAddress + masks_strideWidth;
      const int inBottomLeftMaskAddress = inTopLeftMaskAddress + masks_strideHeight;
      const int inBottomRightMaskAddress = inBottomLeftMaskAddress + masks_strideWidth;

      const int gradMasksTopLeftAddress = gradMasks_strideBatch * b + gradMasks_strideHeight * yInTopLeft + gradMasks_strideWidth * xInTopLeft;
      const int gradMasksTopRightAddress = gradMasksTopLeftAddress + gradMasks_strideWidth;
      const int gradMasksBottomLeftAddress = gradMasksTopLeftAddress + gradMasks_strideHeight;
      const int gradMasksBottomRightAddress = gradMasksBottomLeftAddress + gradMasks_strideWidth;

      const int gradOutputAddress = gradOutput_strideBatch * b + gradOutput_strideHeight * yOut + gradOutput_strideWidth * xOut;

      float topLeftDotProduct = 0;
      float topRightDotProduct = 0;
      float bottomLeftDotProduct = 0;
      float bottomRightDotProduct = 0;

      bool topLeftIsIn = between(xInTopLeft, 0, width-1) && between(yInTopLeft, 0, height-1);
      bool topRightIsIn = between(xInTopLeft+1, 0, width-1) && between(yInTopLeft, 0, height-1);
      bool bottomLeftIsIn = between(xInTopLeft, 0, width-1) && between(yInTopLeft+1, 0, height-1);
      bool bottomRightIsIn = between(xInTopLeft+1, 0, width-1) && between(yInTopLeft+1, 0, height-1);

      float v = 0;
      float inTopLeft=0;
      float inTopRight=0;
      float inBottomLeft=0;
      float inBottomRight=0;

      float c = 0;

      float m = 0;
      float inTopLeftMask=0;
      float inTopRightMask=0;
      float inBottomLeftMask=0;
      float inBottomRightMask=0;

      if(topLeftIsIn) inTopLeftMask = masks_data[inTopLeftMaskAddress];
      if(topRightIsIn) inTopRightMask = masks_data[inTopRightMaskAddress];
      if(bottomLeftIsIn) inBottomLeftMask = masks_data[inBottomLeftMaskAddress];
      if(bottomRightIsIn) inBottomRightMask = masks_data[inBottomRightMaskAddress];


      m = xWeightTopLeft * yWeightTopLeft * inTopLeftMask
        + (1 - xWeightTopLeft) * yWeightTopLeft * inTopRightMask
        + xWeightTopLeft * (1 - yWeightTopLeft) * inBottomLeftMask
        + (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * inBottomRightMask;

      /*
         In that loop we accumulate
         - gradients into the gradInputImages array with atomic adds
         - we compute the dot product that we need for the grid gradient
      */

      for(int t=threadIdx.x; t<inputImages_channels; t+= blockDim.x)
      {

         float gradOutValue = gradOutput_data[gradOutputAddress + t];
         float gradOutValue_fg = gradOutValue * m;
         float gradOutValue_bg = gradOutValue * (1 - m);
         // bool between(int value, int lowerBound, int upperBound)
         if(topLeftIsIn)
         {
            float inTopLeft = inputImages_data[inTopLeftAddress + t];
            topLeftDotProduct += inTopLeft * gradOutValue_fg;
            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesTopLeftAddress + t], xWeightTopLeft * yWeightTopLeft * gradOutValue_fg);
         }

         if(topRightIsIn)
         {
            float inTopRight = inputImages_data[inTopRightAddress + t];
            topRightDotProduct += inTopRight * gradOutValue_fg;
            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesTopRightAddress + t], (1 - xWeightTopLeft) * yWeightTopLeft * gradOutValue_fg);
         }

         if(bottomLeftIsIn)
         {
            float inBottomLeft = inputImages_data[inBottomLeftAddress + t];
            bottomLeftDotProduct += inBottomLeft * gradOutValue_fg;
            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesBottomLeftAddress + t], xWeightTopLeft * (1 - yWeightTopLeft) * gradOutValue_fg);
         }

         if(bottomRightIsIn)
         {
            float inBottomRight = inputImages_data[inBottomRightAddress + t];
            bottomRightDotProduct += inBottomRight * gradOutValue_fg;
            if(!onlyGrid) atomicAdd(&gradInputImages_data[gradInputImagesBottomRightAddress + t], (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * gradOutValue_fg);
         }

         // jw2yang: copy the gradients outside the object region to canvas, and inside region
         if (!topLeftIsIn && !topRightIsIn && !bottomLeftIsIn && !bottomRightIsIn) {
            gradCanvas_data[gradOutputAddress + t] = gradOutValue;
         }
         else {
            gradCanvas_data[gradOutputAddress + t] = gradOutValue_bg;
         }

         // jw2yang: compute the gradient mask value
         if(topLeftIsIn) inTopLeft = inputImages_data[inTopLeftAddress + t];
         if(topRightIsIn) inTopRight = inputImages_data[inTopRightAddress + t];
         if(bottomLeftIsIn) inBottomLeft = inputImages_data[inBottomLeftAddress + t];
         if(bottomRightIsIn) inBottomRight = inputImages_data[inBottomRightAddress + t];
         v = xWeightTopLeft * yWeightTopLeft * inTopLeft
           + (1 - xWeightTopLeft) * yWeightTopLeft * inTopRight
           + xWeightTopLeft * (1 - yWeightTopLeft) * inBottomLeft
           + (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * inBottomRight;

         c = canvas_data[gradOutputAddress + t];

         float gradMaskValue = gradOutValue * (v - c);


         // update gradient on mask map
         if(topLeftIsIn)
         {
            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksTopLeftAddress], xWeightTopLeft * yWeightTopLeft * gradMaskValue);
         }

         if(topRightIsIn)
         {
            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksTopRightAddress], (1 - xWeightTopLeft) * yWeightTopLeft * gradMaskValue);
         }

         if(bottomLeftIsIn)
         {
            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksBottomLeftAddress], xWeightTopLeft * (1 - yWeightTopLeft) * gradMaskValue);
         }

         if(bottomRightIsIn)
         {
            if(!onlyGrid) atomicAdd(&gradMasks_data[gradMasksBottomRightAddress], (1 - xWeightTopLeft) * (1 - yWeightTopLeft) * gradMaskValue);
         }
      }
      /*
         Here we reduce the dot product and compute the grid gradient before writing it.
      */

      /* could do shuffles and use no shmem at all but cuda arch is 2.0 */
      __shared__ volatile float __shmem[16][32];
      __shmem[threadIdx.y][threadIdx.x] = topLeftDotProduct;
      sumReduceShMem(__shmem[threadIdx.y]);
      topLeftDotProduct = __shmem[threadIdx.y][0];

      __shmem[threadIdx.y][threadIdx.x] = topRightDotProduct;
      sumReduceShMem(__shmem[threadIdx.y]);
      topRightDotProduct = __shmem[threadIdx.y][0];

      __shmem[threadIdx.y][threadIdx.x] = bottomLeftDotProduct;
      sumReduceShMem(__shmem[threadIdx.y]);
      bottomLeftDotProduct = __shmem[threadIdx.y][0];

      __shmem[threadIdx.y][threadIdx.x] = bottomRightDotProduct;
      sumReduceShMem(__shmem[threadIdx.y]);
      bottomRightDotProduct = __shmem[threadIdx.y][0];

      yf = - xWeightTopLeft * topLeftDotProduct + xWeightTopLeft * bottomLeftDotProduct - (1-xWeightTopLeft) * topRightDotProduct + (1-xWeightTopLeft) * bottomRightDotProduct;
      xf = - yWeightTopLeft * topLeftDotProduct + yWeightTopLeft * topRightDotProduct - (1-yWeightTopLeft) * bottomLeftDotProduct + (1-yWeightTopLeft) * bottomRightDotProduct;

      if(threadIdx.x==0)
      {
         gridData[threadIdx.y*2] = yf * (inputImages_height-1) / 2;
         gridData[threadIdx.y*2+1] = xf * (inputImages_width-1) / 2;
      }
   }// must put a big if condition in order not to hang at __syncthreads()...
   __syncthreads();

   if(threadIdx.y==0 && withinGridBounds)
       gradGrids_data[b*gradGrids_strideBatch + yOut*gradGrids_strideHeight + xOut*gradGrids_strideWidth + threadIdx.x] = gridData[threadIdx.x];
}





static int cunn_SubSamplerBHWD_updateGradInput(lua_State *L)
{

   float xcoord = (x + 1) * (width - 1) / 2;
   point = floor(xcoord);
   weight = 1 - (xcoord - point);
}

__device__ void getBottomRight(float x, int width, int& point, float& weight)
{
   /* for interpolation :
      stores in point and weight :
      - the x-coordinate of the pixel on the left (or y-coordinate of the upper pixel)
      - the weight for interpolating
   */

   float xcoord = (x + 1) * (width - 1) / 2;
   point = ceil(xcoord);
   weight = 1 - (point - xcoord);
}

__device__ bool between(int value, int lowerBound, int upperBound)
