```python
import torch
import torch.nn as nn
import torch.nn.functional as F
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>
#include <cstdint>

constexpr int THREADS = 256;
constexpr int TILE = 8192;

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

__device__ __forceinline__ Moments combine(Moments a, Moments b) {
    if (a.count == 0) return b;
    if (b.count == 0) return a;
    const int count = a.count + b.count;
    const float delta = b.mean - a.mean;
    const float fraction = float(b.count) / float(count);
    a.mean += delta * fraction;
    a.m2 += b.m2 + delta * delta * float(a.count) * fraction;
    a.count = count;
    return a;
}

__device__ __forceinline__ Moments warp_reduce(Moments value) {
    for (int offset = 16; offset > 0; offset >>= 1) {
        Moments other;
        other.mean = __shfl_down_sync(0xffffffff, value.mean, offset);
        other.m2 = __shfl_down_sync(0xffffffff, value.m2, offset);
        other.count = __shfl_down_sync(0xffffffff, value.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            value = combine(value, other);
    }
    return value;
}

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

    value = warp_reduce(value);
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

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

    value = Moments{0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8)
            value = Moments{means[lane], m2s[lane], counts[lane]};
        value = warp_reduce(value);
    }
    return value;
}

__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial,
    int channels,
    int64_t spatial,
    int64_t elements_per_channel,
    int partitions)
{
    const int channel = blockIdx.y;
    const int partition = blockIdx.x;
    const int64_t begin = int64_t(partition) * TILE;
    const int64_t end =
        begin + TILE < elements_per_channel
            ? begin + TILE
            : elements_per_channel;

    Moments value{0.0f, 0.0f, 0};
    for (int64_t i = begin + threadIdx.x; i < end; i += THREADS) {
        const int64_t batch = i / spatial;
        const int64_t position = i - batch * spatial;
        const float v = x[(batch * channels + channel) * spatial + position];

        ++value.count;
        const float delta = v - value.mean;
        value.mean += delta / float(value.count);
        value.m2 += delta * (v - value.mean);
    }

    value = block_reduce(value);
    if (threadIdx.x == 0) {
        const int64_t index = int64_t(channel) * partitions + partition;
        partial[index * 2] = value.mean;
        partial[index * 2 + 1] = value.m2;
    }
}

__global__ void finish_moments_kernel(
    const float* __restrict__ partial,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ coefficients,
    int64_t elements_per_channel,
    int partitions,
    float momentum,
    float eps)
{
    const int channel = blockIdx.x;
    Moments value{0.0f, 0.0f, 0};

    for (int p = threadIdx.x; p < partitions; p += THREADS) {
        const int64_t remaining = elements_per_channel - int64_t(p) * TILE;
        const int count = remaining < TILE ? int(remaining) : TILE;
        const int64_t index = (int64_t(channel) * partitions + p) * 2;
        value = combine(value, Moments{partial[index], partial[index + 1], count});
    }

    value = block_reduce(value);
    if (threadIdx.x == 0) {
        const float variance = fmaxf(value.m2 / float(elements_per_channel), 0.0f);
        const float invstd = rsqrtf(variance + eps);

        coefficients[channel * 3] = value.mean;
        coefficients[channel * 3 + 1] = invstd * weight[channel];
        coefficients[channel * 3 + 2] = bias[channel];

        running_mean[channel] =
            (1.0f - momentum) * running_mean[channel] + momentum * value.mean;
        running_var[channel] =
            (1.0f - momentum) * running_var[channel] +
            momentum * (value.m2 / float(elements_per_channel - 1));
    }
}

__global__ void eval_coefficients_kernel(
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ running_mean,
    const float* __restrict__ running_var,
    float* __restrict__ coefficients,
    int channels,
    float eps)
{
    const int channel = blockIdx.x * blockDim.x + threadIdx.x;
    if (channel < channels) {
        coefficients[channel * 3] = running_mean[channel];
        coefficients[channel * 3 + 1] =
            weight[channel] * rsqrtf(running_var[channel] + eps);
        coefficients[channel * 3 + 2] = bias[channel];
    }
}

__global__ void pool_normalize_kernel(
    const float* __restrict__ x,
    const float* __restrict__ coefficients,
    float* __restrict__ output,
    int channels,
    int depth,
    int height,
    int width,
    int out_depth,
    int out_height,
    int out_width,
    int64_t output_elements)
{
    const int64_t index = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    if (index >= output_elements) return;

    int64_t remaining = index;
    const int ow = remaining % out_width;
    remaining /= out_width;
    const int oh = remaining % out_height;
    remaining /= out_height;
    const int od = remaining % out_depth;
    remaining /= out_depth;
    const int channel = remaining % channels;

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

    float sum = 0.0f;
    #pragma unroll
    for (int z = 0; z < 4; ++z) {
        #pragma unroll
        for (int y = 0; y < 4; ++y) {
            const int64_t offset = base + int64_t(z) * height * width + y * width;
            const float row =
                (x[offset] + x[offset + 1]) +
                (x[offset + 2] + x[offset + 3]);
            sum += row;
        }
    }

    const float mean = coefficients[channel * 3];
    const float scale = coefficients[channel * 3 + 1];
    const float bias = coefficients[channel * 3 + 2];
    output[index] = (sum * (1.0f / 64.0f) - mean) * scale + bias;
}

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)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == torch::kFloat32,
                "Expected a CUDA float32 input");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");

    const c10::cuda::CUDAGuard device_guard(x.device());
    const int batches = x.size(0);
    const int channels = x.size(1);
    const int depth = x.size(2);
    const int height = x.size(3);
    const int width = x.size(4);
    TORCH_CHECK(depth >= 4 && height >= 4 && width >= 4,
                "Spatial dimensions must be at least four");

    const int out_depth = depth / 4;
    const int out_height = height / 4;
    const int out_width = width / 4;
    const int64_t spatial = int64_t(depth) * height * width;
    const int64_t elements_per_channel = int64_t(batches) * spatial;

    auto output = torch::empty(
        {batches, channels, out_depth, out_height, out_width}, x.options());
    if (batches == 0) return output;

    auto coefficients = torch::empty({channels, 3}, x.options());
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        TORCH_CHECK(elements_per_channel > 1,
                    "Batch normalization requires more than one value per channel");
        const int partitions = (elements_per_channel + TILE - 1) / TILE;
        auto partial = torch::empty({channels, partitions, 2}, x.options());

        partial_moments_kernel<<<dim3(partitions, channels), THREADS, 0, stream>>>(
            x.data_ptr<float>(), partial.data_ptr<float>(), channels,
            spatial, elements_per_channel, partitions);
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        finish_moments_kernel<<<channels, THREADS, 0, stream>>>(
            partial.data_ptr<float>(),
            weight.data_ptr<float>(), bias.data_ptr<float>(),
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            coefficients.data_ptr<float>(), elements_per_channel, partitions,
            float(momentum), float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    } else {
        eval_coefficients_kernel<<<(channels + THREADS - 1) / THREADS,
                                   THREADS, 0, stream>>>(
            weight.data_ptr<float>(), bias.data_ptr<float>(),
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            coefficients.data_ptr<float>(), channels, float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    const int64_t output_elements = output.numel();
    pool_normalize_kernel<<<(output_elements + THREADS - 1) / THREADS,
                            THREADS, 0, stream>>>(
        x.data_ptr<float>(), coefficients.data_ptr<float>(),
        output.data_ptr<float>(), channels, depth, height, width,
        out_depth, out_height, out_width, output_elements);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_extension = None


def _get_extension():
    global _extension
    if _extension is None:
        _extension = 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,
        )
    return _extension


class _BatchNormPool(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, running_mean, running_var,
                training, momentum, eps):
        ctx.training = training
        ctx.eps = eps

        if any(ctx.needs_input_grad[:3]):
            if training:
                ctx.save_for_backward(x, weight, bias)
            else:
                ctx.save_for_backward(
                    x, weight, bias, running_mean.clone(), running_var.clone()
                )

        return _get_extension().bn_pool_cuda(
            x, weight, bias, running_mean, running_var,
            training, momentum, eps
        )

    @staticmethod
    def backward(ctx, grad_output):
        saved = ctx.saved_tensors
        required = ctx.needs_input_grad[:3]

        with torch.enable_grad():
            x = saved[0].detach().requires_grad_(required[0])
            weight = saved[1].detach().requires_grad_(required[1])
            bias = saved[2].detach().requires_grad_(required[2])
            running_mean = None if ctx.training else saved[3]
            running_var = None if ctx.training else saved[4]

            normalized = F.batch_norm(
                x, running_mean, running_var, weight, bias,
                training=ctx.training, momentum=0.0, eps=ctx.eps
            )
            output = F.avg_pool3d(F.avg_pool3d(normalized, 2), 2)
            targets = [
                tensor for tensor, needed in zip((x, weight, bias), required)
                if needed
            ]
            computed = torch.autograd.grad(output, targets, grad_output)

        gradients = iter(computed)
        dx, dw, db = [
            next(gradients) if needed else None for needed in required
        ]
        return dx, dw, db, None, None, None, None, None


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

        supported = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and bn.affine
            and bn.track_running_stats
            and all(
                tensor is not None
                and tensor.dtype == torch.float32
                and tensor.is_contiguous()
                and tensor.device == x.device
                for tensor in (
                    bn.weight, bn.bias, bn.running_mean, bn.running_var
                )
            )
            and min(x.shape[2:]) >= 4
        )
        if not supported:
            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 _BatchNormPool.apply(
            x, bn.weight, bn.bias, bn.running_mean, bn.running_var,
            bn.training, momentum, bn.eps
        )
```