```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,
    torch::Tensor batches,
    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 <cstdint>

constexpr int THREADS = 256;
constexpr int ITEMS = 32;
constexpr int CHUNK = THREADS * ITEMS;

__device__ __forceinline__ void reduce_pair(float& a, float& b) {
    __shared__ float sa[8];
    __shared__ float sb[8];

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

    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        a += __shfl_down_sync(0xffffffff, a, offset);
        b += __shfl_down_sync(0xffffffff, b, offset);
    }

    if (lane == 0) {
        sa[warp] = a;
        sb[warp] = b;
    }
    __syncthreads();

    if (warp == 0) {
        a = lane < 8 ? sa[lane] : 0.0f;
        b = lane < 8 ? sb[lane] : 0.0f;

        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1) {
            a += __shfl_down_sync(0xffffffff, a, offset);
            b += __shfl_down_sync(0xffffffff, b, offset);
        }
    }
}

__global__ void moments_kernel(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int64_t* batches,
    int64_t spatial,
    int channels,
    int batch_size,
    int chunks)
{
    const int chunk = blockIdx.x;
    const int n = blockIdx.y;
    const int c = blockIdx.z;

    if (chunk == 0 && n == 0 && c == 0 && threadIdx.x == 0) {
        batches[0] += 1;
    }

    const int64_t base = (int64_t(n) * channels + c) * spatial;
    const int64_t start = int64_t(chunk) * CHUNK + threadIdx.x;

    float sum = 0.0f;
    float sum_sq = 0.0f;

    #pragma unroll
    for (int i = 0; i < ITEMS; ++i) {
        const int64_t pos = start + i * THREADS;
        if (pos < spatial) {
            const float v = x[base + pos];
            sum += v;
            sum_sq = fmaf(v, v, sum_sq);
        }
    }

    reduce_pair(sum, sum_sq);

    if (threadIdx.x == 0) {
        const int64_t index =
            (int64_t(c) * batch_size + n) * chunks + chunk;
        partial[index] = make_float2(sum, sum_sq);
    }
}

__global__ void finalize_kernel(
    const float2* __restrict__ partial,
    float* __restrict__ means,
    float* __restrict__ invstds,
    float* running_mean,
    float* running_var,
    const int64_t* batches,
    int partial_count,
    int64_t sample_count,
    float momentum,
    float eps)
{
    const int c = blockIdx.x;
    float sum = 0.0f;
    float sum_sq = 0.0f;

    for (int i = threadIdx.x; i < partial_count; i += THREADS) {
        const float2 p = partial[int64_t(c) * partial_count + i];
        sum += p.x;
        sum_sq += p.y;
    }

    reduce_pair(sum, sum_sq);

    if (threadIdx.x == 0) {
        const double count = double(sample_count);
        const double mean_d = double(sum) / count;
        const float mean = float(mean_d);
        const float var = float(fmax(
            double(sum_sq) / count - mean_d * mean_d, 0.0));

        means[c] = mean;
        invstds[c] = rsqrtf(var + eps);

        const float factor = momentum < 0.0f
            ? float(1.0 / double(batches[0]))
            : momentum;

        const float unbiased_var =
            float(double(var) * count / (count - 1.0));

        running_mean[c] =
            (1.0f - factor) * running_mean[c] + factor * mean;
        running_var[c] =
            (1.0f - factor) * running_var[c] + factor * unbiased_var;
    }
}

template <bool Training>
__global__ void pool_normalize_kernel(
    const float* __restrict__ x,
    const float* __restrict__ means,
    const float* __restrict__ variance_or_invstd,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ out,
    int64_t output_count,
    int channels,
    int depth,
    int height,
    int width,
    int out_depth,
    int out_height,
    int out_width,
    float eps)
{
    // Four neighboring threads cooperatively read one 4x4x4 cube.
    const int lane = threadIdx.x & 3;
    const int64_t index =
        int64_t(blockIdx.x) * (THREADS / 4) + (threadIdx.x >> 2);

    float sum = 0.0f;
    int c = 0;

    if (index < output_count) {
        int64_t q = index;
        const int ow = q % out_width;
        q /= out_width;
        const int oh = q % out_height;
        q /= out_height;
        const int od = q % out_depth;
        const int64_t nc = q / out_depth;
        c = nc % channels;

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

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

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

    if (index < output_count && lane == 0) {
        const float invstd = Training
            ? variance_or_invstd[c]
            : rsqrtf(variance_or_invstd[c] + eps);

        const float normalized = (sum * (1.0f / 64.0f) - means[c]) * invstd;
        out[index] = normalized * weight[c] + bias[c];
    }
}

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

    const int n = x.size(0);
    const int c = x.size(1);
    const int d = x.size(2);
    const int h = x.size(3);
    const int w = x.size(4);

    const int od = d / 4;
    const int oh = h / 4;
    const int ow = w / 4;
    const int64_t spatial = int64_t(d) * h * w;

    auto out = torch::empty({n, c, od, oh, ow}, x.options());
    const int64_t output_count = out.numel();
    const int pool_blocks = (output_count + THREADS / 4 - 1) / (THREADS / 4);

    if (training) {
        const int chunks = (spatial + CHUNK - 1) / CHUNK;
        const int partial_count = n * chunks;

        auto partial = torch::empty({c, partial_count, 2}, x.options());
        auto stats = torch::empty({2, c}, x.options());
        float* means = stats.data_ptr<float>();
        float* invstds = means + c;

        moments_kernel<<<dim3(chunks, n, c), THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            batches.data_ptr<int64_t>(),
            spatial, c, n, chunks);
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        finalize_kernel<<<c, THREADS, 0, stream>>>(
            reinterpret_cast<const float2*>(partial.data_ptr<float>()),
            means,
            invstds,
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            batches.data_ptr<int64_t>(),
            partial_count,
            int64_t(n) * spatial,
            float(momentum),
            float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        pool_normalize_kernel<true><<<pool_blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            means,
            invstds,
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            out.data_ptr<float>(),
            output_count,
            c, d, h, w, od, oh, ow,
            float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    } else {
        pool_normalize_kernel<false><<<pool_blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            out.data_ptr<float>(),
            output_count,
            c, d, h, w, od, oh, ow,
            float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    return out;
}
"""

_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 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

        use_custom = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 5
            and x.shape[0] > 0
            and min(x.shape[2:]) >= 4
            and x.shape[0] <= 65535
            and x.shape[1] <= 65535
            and bn.affine
            and bn.track_running_stats
            and bn.running_mean is not None
            and bn.running_var is not None
            and bn.num_batches_tracked is not None
            and bn.weight.dtype == torch.float32
            and bn.bias.dtype == torch.float32
            and bn.running_mean.dtype == torch.float32
            and bn.running_var.dtype == torch.float32
            and bn.num_batches_tracked.is_cuda
        )

        if not use_custom:
            return self.avg_pool2(self.avg_pool1(bn(x)))

        return _fused_ops.fused_bn_pool_cuda(
            x.contiguous(),
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            bn.num_batches_tracked,
            bn.training,
            -1.0 if bn.momentum is None else bn.momentum,
            bn.eps,
        )
```