```python
import torch
import torch.nn as nn
from torch.utils.cpp_extension import load_inline


_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    double eps);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <algorithm>

struct Moments {
    float mean;
    float m2;
    int count;
};

__device__ __forceinline__ Moments merge_moments(Moments a, Moments b) {
    if (b.count == 0) return a;
    if (a.count == 0) return b;

    int count = a.count + b.count;
    float delta = b.mean - a.mean;
    float ratio = float(b.count) / float(count);

    a.mean += delta * ratio;
    a.m2 += b.m2 + delta * delta * float(a.count) * ratio;
    a.count = count;
    return a;
}

__device__ __forceinline__ Moments warp_reduce(Moments a) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Moments b;
        b.mean = __shfl_down_sync(0xffffffff, a.mean, offset);
        b.m2 = __shfl_down_sync(0xffffffff, a.m2, offset);
        b.count = __shfl_down_sync(0xffffffff, a.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            a = merge_moments(a, b);
    }
    return a;
}

__device__ __forceinline__ Moments block_reduce(Moments a) {
    __shared__ float means[8];
    __shared__ float m2s[8];
    __shared__ int counts[8];

    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;

    a = warp_reduce(a);
    if (lane == 0) {
        means[warp] = a.mean;
        m2s[warp] = a.m2;
        counts[warp] = a.count;
    }
    __syncthreads();

    Moments result = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8) {
            result.mean = means[lane];
            result.m2 = m2s[lane];
            result.count = counts[lane];
        }
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_mean,
    float* __restrict__ partial_m2,
    int* __restrict__ partial_count,
    int channels,
    int spatial,
    int channel_elements,
    int chunks)
{
    int channel = blockIdx.y;
    Moments a = {0.0f, 0.0f, 0};

    for (int k = blockIdx.x * blockDim.x + threadIdx.x;
         k < channel_elements;
         k += chunks * blockDim.x) {
        int batch = k / spatial;
        int position = k - batch * spatial;
        int64_t index =
            (int64_t(batch) * channels + channel) * spatial + position;

        float value = x[index];
        ++a.count;
        float delta = value - a.mean;
        a.mean += delta / float(a.count);
        a.m2 += delta * (value - a.mean);
    }

    a = block_reduce(a);
    if (threadIdx.x == 0) {
        int index = channel * chunks + blockIdx.x;
        partial_mean[index] = a.mean;
        partial_m2[index] = a.m2;
        partial_count[index] = a.count;
    }
}

__global__ void finish_moments_kernel(
    const float* __restrict__ partial_mean,
    const float* __restrict__ partial_m2,
    const int* __restrict__ partial_count,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ stats,
    int channels,
    int chunks,
    float momentum,
    float eps)
{
    int channel = blockIdx.x;
    Moments a = {0.0f, 0.0f, 0};

    for (int k = threadIdx.x; k < chunks; k += blockDim.x) {
        int index = channel * chunks + k;
        Moments b = {
            partial_mean[index],
            partial_m2[index],
            partial_count[index]
        };
        a = merge_moments(a, b);
    }

    a = block_reduce(a);
    if (threadIdx.x == 0) {
        float variance = fmaxf(a.m2 / float(a.count), 0.0f);
        stats[channel] = a.mean;
        stats[channels + channel] = rsqrtf(variance + eps);

        running_mean[channel] =
            (1.0f - momentum) * running_mean[channel] + momentum * a.mean;
        running_var[channel] =
            (1.0f - momentum) * running_var[channel]
            + momentum * variance * (float(a.count) / float(a.count - 1));
    }
}

__global__ void evaluation_stats_kernel(
    const float* __restrict__ running_mean,
    const float* __restrict__ running_var,
    float* __restrict__ stats,
    int channels,
    float eps)
{
    int channel = blockIdx.x * blockDim.x + threadIdx.x;
    if (channel < channels) {
        stats[channel] = running_mean[channel];
        stats[channels + channel] = rsqrtf(running_var[channel] + eps);
    }
}

// Four neighboring lanes cooperate on each 4x4x4 pooling region.
// Two successive non-overlapping average pools are a single average
// over this region. Batch normalization commutes with that average.
__global__ void normalized_pool_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ stats,
    float* __restrict__ output,
    int channels,
    int depth,
    int height,
    int width,
    int output_depth,
    int output_height,
    int output_width,
    int64_t output_elements)
{
    int64_t thread = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    int64_t output_index = thread >> 2;
    int lane = threadIdx.x & 3;

    float sum = 0.0f;
    int channel = 0;

    if (output_index < output_elements) {
        int64_t remaining = output_index;
        int ow = remaining % output_width;
        remaining /= output_width;
        int oh = remaining % output_height;
        remaining /= output_height;
        int od = remaining % output_depth;
        remaining /= output_depth;
        channel = remaining % channels;

        int64_t base =
            ((remaining * depth + od * 4) * height + oh * 4) * width
            + ow * 4 + lane;

        #pragma unroll
        for (int d = 0; d < 4; ++d) {
            #pragma unroll
            for (int h = 0; h < 4; ++h) {
                sum += x[base + int64_t(d) * height * width + h * width];
            }
        }
    }

    sum += __shfl_down_sync(0xffffffff, sum, 2, 4);
    sum += __shfl_down_sync(0xffffffff, sum, 1, 4);

    if (lane == 0 && output_index < output_elements) {
        float average = sum * (1.0f / 64.0f);
        float normalized =
            (average - stats[channel]) * stats[channels + channel];
        output[output_index] = normalized * weight[channel] + bias[channel];
    }
}

torch::Tensor bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    double eps)
{
    c10::cuda::CUDAGuard guard(x.device());
    auto stream = at::cuda::getCurrentCUDAStream();

    int batches = x.size(0);
    int channels = x.size(1);
    int depth = x.size(2);
    int height = x.size(3);
    int width = x.size(4);
    int spatial = depth * height * width;
    int channel_elements = batches * spatial;

    int od = depth / 4;
    int oh = height / 4;
    int ow = width / 4;

    auto output = torch::empty({batches, channels, od, oh, ow}, x.options());
    auto stats = torch::empty({2, channels}, x.options());

    if (training) {
        int chunks = std::min(1024, (channel_elements + 4095) / 4096);
        chunks = std::max(chunks, 1);

        auto partial_mean = torch::empty({channels, chunks}, x.options());
        auto partial_m2 = torch::empty({channels, chunks}, x.options());
        auto partial_count = torch::empty(
            {channels, chunks}, x.options().dtype(torch::kInt32));

        partial_moments_kernel<<<dim3(chunks, channels), 256, 0, stream>>>(
            x.data_ptr<float>(),
            partial_mean.data_ptr<float>(),
            partial_m2.data_ptr<float>(),
            partial_count.data_ptr<int>(),
            channels, spatial, channel_elements, chunks);
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        finish_moments_kernel<<<channels, 256, 0, stream>>>(
            partial_mean.data_ptr<float>(),
            partial_m2.data_ptr<float>(),
            partial_count.data_ptr<int>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            stats.data_ptr<float>(),
            channels, chunks, float(momentum), float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    } else {
        evaluation_stats_kernel<<<(channels + 255) / 256, 256, 0, stream>>>(
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            stats.data_ptr<float>(),
            channels, float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    int64_t elements = output.numel();
    if (elements != 0) {
        normalized_pool_kernel<<<(elements + 63) / 64, 256, 0, stream>>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            stats.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, depth, height, width, od, oh, ow, elements);
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    return output;
}
"""

_fused_bn_pool = load_inline(
    name="fused_bn_double_avgpool3d_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["bn_pool_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_channels, out_channels, kernel_size, stride, padding, bias_shape
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
        )
        self.batch_norm = nn.BatchNorm3d(out_channels)
        self.avg_pool1 = nn.AvgPool3d(kernel_size=2)
        self.avg_pool2 = nn.AvgPool3d(kernel_size=2)

    def forward(self, x):
        x = self.conv_transpose(x)
        bn = self.batch_norm

        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or not x.is_contiguous()
            or min(x.shape[2:]) < 4
            or not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean.dtype != torch.float32
        ):
            return self.avg_pool2(self.avg_pool1(bn(x)))

        momentum = bn.momentum if bn.momentum is not None else 0.0
        if bn.training:
            bn.num_batches_tracked.add_(1)
            if bn.momentum is None:
                momentum = 1.0 / float(bn.num_batches_tracked.item())

        return _fused_bn_pool.bn_pool_cuda(
            x,
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            bn.training,
            momentum,
            bn.eps,
        )
```