```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 fused_bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double eps,
    double momentum,
    bool training);
"""

_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 <cmath>
#include <cstdint>

constexpr int THREADS = 256;
constexpr int CHUNK = 4096;

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

__device__ __forceinline__ Moments merge_moments(
    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 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)
{
    const int lane = threadIdx.x & 31;
    #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 (lane + offset < 32)
            a = merge_moments(a, b);
    }
    return a;
}

__device__ __forceinline__ Moments block_reduce(Moments a)
{
    __shared__ Moments warp_results[THREADS / 32];

    a = warp_reduce(a);
    if ((threadIdx.x & 31) == 0)
        warp_results[threadIdx.x >> 5] = a;
    __syncthreads();

    Moments result = {0.0f, 0.0f, 0};
    if (threadIdx.x < 32) {
        if (threadIdx.x < THREADS / 32)
            result = warp_results[threadIdx.x];
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int64_t spatial,
    int64_t elements_per_channel,
    int channels,
    int chunks)
{
    const int channel = blockIdx.y;
    const int chunk = blockIdx.x;
    const int64_t start = int64_t(chunk) * CHUNK;

    Moments local = {0.0f, 0.0f, 0};

    #pragma unroll
    for (int k = 0; k < CHUNK / THREADS; ++k) {
        const int64_t j = start + threadIdx.x + k * THREADS;
        if (j < elements_per_channel) {
            const int64_t batch = j / spatial;
            const int64_t position = j - batch * spatial;
            const float value =
                x[(batch * channels + channel) * spatial + position];

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

    Moments total = block_reduce(local);
    if (threadIdx.x == 0)
        partial[int64_t(channel) * chunks + chunk] =
            make_float2(total.mean, total.m2);
}

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

    for (int chunk = threadIdx.x; chunk < chunks; chunk += THREADS) {
        const float2 value = partial[int64_t(channel) * chunks + chunk];
        const int64_t remaining =
            elements_per_channel - int64_t(chunk) * CHUNK;
        const int count = remaining < CHUNK ? int(remaining) : CHUNK;
        Moments part = {value.x, value.y, count};
        local = merge_moments(local, part);
    }

    Moments total = block_reduce(local);
    if (threadIdx.x == 0) {
        const float variance = fmaxf(total.m2 / float(total.count), 0.0f);
        const float scale = weight[channel] * rsqrtf(variance + eps);
        coefficients[channel] =
            make_float2(scale, bias[channel] - total.mean * scale);

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

__global__ void evaluation_coefficients_kernel(
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ running_mean,
    const float* __restrict__ running_var,
    float2* __restrict__ coefficients,
    int channels,
    float eps)
{
    const int channel = blockIdx.x * blockDim.x + threadIdx.x;
    if (channel < channels) {
        const float scale =
            weight[channel] * rsqrtf(running_var[channel] + eps);
        coefficients[channel] =
            make_float2(scale, bias[channel] - running_mean[channel] * scale);
    }
}

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

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

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

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

    const float2 affine = coefficients[channel];
    output[index] = fmaf(sum * (1.0f / 64.0f), affine.x, affine.y);
}

torch::Tensor fused_bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double eps,
    double momentum,
    bool training)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == torch::kFloat32,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");

    const c10::cuda::CUDAGuard device_guard(x.device());
    const auto stream = at::cuda::getCurrentCUDAStream();

    const int64_t batch = 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);
    const int out_depth = depth / 4;
    const int out_height = height / 4;
    const int out_width = width / 4;

    TORCH_CHECK(out_depth > 0 && out_height > 0 && out_width > 0,
                "Input dimensions must be at least four");

    const int64_t spatial = int64_t(depth) * height * width;
    const int64_t elements_per_channel = batch * spatial;

    auto output = torch::empty(
        {batch, channels, out_depth, out_height, out_width}, x.options());
    auto coefficients = torch::empty({channels, 2}, x.options());

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

        partial_moments_kernel<<<dim3(chunks, channels), THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, elements_per_channel, channels, chunks);

        training_coefficients_kernel<<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const float2*>(partial.data_ptr<float>()),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            reinterpret_cast<float2*>(coefficients.data_ptr<float>()),
            elements_per_channel, chunks, float(eps), float(momentum));
    } else {
        evaluation_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>(),
            reinterpret_cast<float2*>(coefficients.data_ptr<float>()),
            channels, float(eps));
    }

    const int64_t output_elements = output.numel();
    if (output_elements > 0) {
        pooled_affine_kernel<<<
            (output_elements + THREADS - 1) / THREADS, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<const float2*>(coefficients.data_ptr<float>()),
            output.data_ptr<float>(),
            output_elements, channels, depth, height, width,
            out_depth, out_height, out_width);
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

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


class _FusedBNPool(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, running_mean, running_var,
                eps, momentum, training):
        output = _fused_ops.fused_bn_pool_cuda(
            x, weight, bias, running_mean, running_var,
            eps, momentum, training,
        )
        ctx.training_mode = training
        ctx.eps = eps
        ctx.save_for_backward(x, weight, running_mean, running_var)
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, running_mean, running_var = ctx.saved_tensors
        axes = (0, 2, 3, 4)
        shape = (1, -1, 1, 1, 1)

        if ctx.training_mode:
            variance, mean = torch.var_mean(x, dim=axes, unbiased=False)
            inv_std = torch.rsqrt(variance + ctx.eps)
        else:
            mean = running_mean
            inv_std = torch.rsqrt(running_var + ctx.eps)

        grad_pooled = grad_output / 64.0
        for dim in (2, 3, 4):
            grad_pooled = grad_pooled.repeat_interleave(4, dim=dim)
        grad_pooled = torch.nn.functional.pad(
            grad_pooled,
            (
                0, x.shape[4] - grad_pooled.shape[4],
                0, x.shape[3] - grad_pooled.shape[3],
                0, x.shape[2] - grad_pooled.shape[2],
            ),
        )

        normalized = (x - mean.view(shape)) * inv_std.view(shape)
        grad_bias = grad_pooled.sum(dim=axes)
        grad_weight = (grad_pooled * normalized).sum(dim=axes)

        if ctx.training_mode:
            count = x.numel() // x.shape[1]
            grad_x = (
                grad_pooled
                - grad_bias.view(shape) / count
                - normalized * (grad_weight.view(shape) / count)
            ) * (weight * inv_std).view(shape)
        else:
            grad_x = grad_pooled * (weight * inv_std).view(shape)

        return grad_x, grad_weight, grad_bias, 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

        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or not bn.affine
            or not bn.track_running_stats
            or any(t.dtype != torch.float32 for t in (
                bn.weight, bn.bias, bn.running_mean, bn.running_var
            ))
        ):
            return self.avg_pool2(self.avg_pool1(bn(x)))

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

        return _FusedBNPool.apply(
            x.contiguous(),
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            bn.eps,
            0.0 if momentum is None else momentum,
            bn.training,
        )
```