```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_pool4_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>

constexpr int THREADS = 256;
constexpr int TILE = 4096;

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 fraction = float(b.count) / float(count);

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

__device__ __forceinline__ Moments warp_reduce(Moments value) {
    int lane = threadIdx.x & 31;
    #pragma unroll
    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 (lane + offset < 32)
            value = merge_moments(value, other);
    }
    return value;
}

__device__ __forceinline__ Moments block_reduce(Moments value) {
    __shared__ Moments warp_values[THREADS / 32];

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

    if (lane == 0)
        warp_values[warp] = value;
    __syncthreads();

    Moments result = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < THREADS / 32)
            result = warp_values[lane];
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_mean,
    float* __restrict__ partial_m2,
    int64_t spatial,
    int64_t samples,
    int channels,
    int parts)
{
    int channel = blockIdx.y;
    int part = blockIdx.x;
    int64_t start = int64_t(part) * TILE + threadIdx.x;

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

    #pragma unroll
    for (int j = 0; j < TILE / THREADS; ++j) {
        int64_t index = start + int64_t(j) * THREADS;
        if (index < samples) {
            int64_t batch = index / spatial;
            int64_t position = index - batch * spatial;
            float v = x[(batch * channels + channel) * spatial + position];

            ++value.count;
            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) {
        int64_t index = int64_t(channel) * parts + part;
        partial_mean[index] = value.mean;
        partial_m2[index] = value.m2;
    }
}

__global__ void finish_moments_kernel(
    const float* __restrict__ partial_mean,
    const float* __restrict__ partial_m2,
    float* __restrict__ means,
    float* __restrict__ invstd,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    int64_t samples,
    int parts,
    float momentum,
    float eps)
{
    int channel = blockIdx.x;
    Moments value = {0.0f, 0.0f, 0};

    for (int part = threadIdx.x; part < parts; part += THREADS) {
        int64_t remaining = samples - int64_t(part) * TILE;
        int count = remaining < TILE ? int(remaining) : TILE;
        int64_t index = int64_t(channel) * parts + part;
        Moments other = {
            partial_mean[index], partial_m2[index], count
        };
        value = merge_moments(value, other);
    }

    value = block_reduce(value);
    if (threadIdx.x == 0) {
        float variance = fmaxf(value.m2 / float(samples), 0.0f);
        means[channel] = value.mean;
        invstd[channel] = rsqrtf(variance + eps);

        running_mean[channel] =
            (1.0f - momentum) * running_mean[channel]
            + momentum * value.mean;

        float unbiased_variance =
            fmaxf(value.m2 / float(samples - 1), 0.0f);
        running_var[channel] =
            (1.0f - momentum) * running_var[channel]
            + momentum * unbiased_variance;
    }
}

template<bool TRAINING>
__global__ void pool_normalize_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ means,
    const float* __restrict__ variance_or_invstd,
    float* __restrict__ output,
    int channels,
    int depth,
    int height,
    int width,
    int out_depth,
    int out_height,
    int out_width,
    int64_t output_count,
    float eps)
{
    int64_t index = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    if (index >= output_count)
        return;

    int64_t remaining = index;
    int ox = remaining % out_width;
    remaining /= out_width;
    int oy = remaining % out_height;
    remaining /= out_height;
    int oz = remaining % out_depth;
    int64_t nc = remaining / out_depth;
    int channel = nc % channels;

    int64_t spatial = int64_t(depth) * height * width;
    int64_t base = nc * spatial
        + int64_t(oz * 4) * height * width
        + int64_t(oy * 4) * width
        + ox * 4;

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

    float scale;
    if (TRAINING)
        scale = variance_or_invstd[channel];
    else
        scale = rsqrtf(variance_or_invstd[channel] + eps);

    float normalized = (sum * (1.0f / 64.0f) - means[channel]) * scale;
    output[index] = normalized * weight[channel] + bias[channel];
}

torch::Tensor bn_pool4_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());
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    int64_t batch = x.size(0);
    int channels = int(x.size(1));
    int depth = int(x.size(2));
    int height = int(x.size(3));
    int width = int(x.size(4));

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

    auto output = torch::empty({batch, channels, od, oh, ow}, x.options());
    if (output.numel() == 0)
        return output;

    torch::Tensor means = running_mean;
    torch::Tensor variance_or_invstd = running_var;

    if (training) {
        int64_t spatial = int64_t(depth) * height * width;
        int64_t samples = batch * spatial;
        TORCH_CHECK(samples > 1, "Batch normalization requires multiple samples");

        int parts = int((samples + TILE - 1) / TILE);
        auto partial_mean = torch::empty({channels, parts}, x.options());
        auto partial_m2 = torch::empty({channels, parts}, x.options());
        means = torch::empty({channels}, x.options());
        variance_or_invstd = torch::empty({channels}, x.options());

        partial_moments_kernel<<<dim3(parts, channels), THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            partial_mean.data_ptr<float>(),
            partial_m2.data_ptr<float>(),
            spatial, samples, channels, parts);
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        finish_moments_kernel<<<channels, THREADS, 0, stream>>>(
            partial_mean.data_ptr<float>(),
            partial_m2.data_ptr<float>(),
            means.data_ptr<float>(),
            variance_or_invstd.data_ptr<float>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            samples, parts, float(momentum), float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    int64_t count = output.numel();
    int blocks = int((count + THREADS - 1) / THREADS);

    if (training) {
        pool_normalize_kernel<true><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(), weight.data_ptr<float>(), bias.data_ptr<float>(),
            means.data_ptr<float>(), variance_or_invstd.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, depth, height, width, od, oh, ow, count, float(eps));
    } else {
        pool_normalize_kernel<false><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(), weight.data_ptr<float>(), bias.data_ptr<float>(),
            means.data_ptr<float>(), variance_or_invstd.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, depth, height, width, od, oh, ow, count, float(eps));
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="convtranspose3d_bn_pool4_fused_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["bn_pool4_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 bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean.dtype != torch.float32
            or min(x.shape[2:]) < 4
        ):
            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 / float(bn.num_batches_tracked.item())

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