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


_cpp_source = r"""
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.h>
#include <cuda_runtime.h>
#include <stdint.h>
#include <cmath>

constexpr int THREADS = 256;
constexpr int CHUNK = 8192;

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

__device__ __forceinline__ Welford merge_stats(Welford a, Welford b) {
    if (b.count == 0) return a;
    if (a.count == 0) return b;
    float delta = b.mean - a.mean;
    int count = a.count + b.count;
    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__ Welford warp_reduce(Welford v) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Welford other;
        other.mean = __shfl_down_sync(0xffffffff, v.mean, offset);
        other.m2 = __shfl_down_sync(0xffffffff, v.m2, offset);
        other.count = __shfl_down_sync(0xffffffff, v.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            v = merge_stats(v, other);
    }
    return v;
}

__device__ __forceinline__ Welford block_reduce(Welford v) {
    __shared__ Welford shared[8];
    v = warp_reduce(v);
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    if (lane == 0) shared[warp] = v;
    __syncthreads();

    Welford result = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8) result = shared[lane];
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_stats_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_mean,
    float* __restrict__ partial_m2,
    int channels,
    int64_t spatial,
    int64_t samples,
    int parts)
{
    int part = blockIdx.x;
    int c = blockIdx.y;
    int64_t begin = int64_t(part) * CHUNK;

    Welford v = {0.0f, 0.0f, 0};

    #pragma unroll
    for (int k = 0; k < CHUNK / THREADS; ++k) {
        int64_t index = begin + threadIdx.x + k * THREADS;
        if (index < samples) {
            int64_t n = index / spatial;
            int64_t s = index - n * spatial;
            float value = x[(n * channels + c) * spatial + s];
            ++v.count;
            float delta = value - v.mean;
            v.mean += delta / float(v.count);
            v.m2 += delta * (value - v.mean);
        }
    }

    v = block_reduce(v);
    if (threadIdx.x == 0) {
        int64_t index = int64_t(c) * parts + part;
        partial_mean[index] = v.mean;
        partial_m2[index] = v.m2;
    }
}

__global__ void finalize_stats_kernel(
    const float* __restrict__ partial_mean,
    const float* __restrict__ partial_m2,
    const float* __restrict__ weight,
    float* running_mean,
    float* running_var,
    float* __restrict__ mean,
    float* __restrict__ scale,
    int64_t samples,
    int parts,
    float momentum,
    float eps)
{
    int c = blockIdx.x;
    Welford v = {0.0f, 0.0f, 0};

    for (int p = threadIdx.x; p < parts; p += THREADS) {
        int64_t remaining = samples - int64_t(p) * CHUNK;
        Welford other;
        other.mean = partial_mean[int64_t(c) * parts + p];
        other.m2 = partial_m2[int64_t(c) * parts + p];
        other.count = remaining < CHUNK ? int(remaining) : CHUNK;
        v = merge_stats(v, other);
    }

    v = block_reduce(v);
    if (threadIdx.x == 0) {
        float variance = fmaxf(v.m2 / float(samples), 0.0f);
        mean[c] = v.mean;
        scale[c] = weight[c] * rsqrtf(variance + eps);

        if (running_mean != nullptr) {
            running_mean[c] =
                (1.0f - momentum) * running_mean[c] + momentum * v.mean;
            float unbiased_variance = v.m2 / float(samples - 1);
            running_var[c] =
                (1.0f - momentum) * running_var[c]
                + momentum * unbiased_variance;
        }
    }
}

__global__ void eval_stats_kernel(
    const float* __restrict__ running_mean,
    const float* __restrict__ running_var,
    const float* __restrict__ weight,
    float* __restrict__ mean,
    float* __restrict__ scale,
    int channels,
    float eps)
{
    int c = blockIdx.x * blockDim.x + threadIdx.x;
    if (c < channels) {
        mean[c] = running_mean[c];
        scale[c] = weight[c] * rsqrtf(running_var[c] + eps);
    }
}

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

    int64_t q = index;
    int ow = q % out_width;
    q /= out_width;
    int oh = q % out_height;
    q /= out_height;
    int od = q % out_depth;
    q /= out_depth;
    int c = q % channels;

    int64_t spatial = int64_t(depth) * height * width;
    int64_t base = q * spatial
        + int64_t(od * 4) * height * width
        + int64_t(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) {
            int64_t row = base + int64_t(dz) * height * width + dy * width;
            float row_sum = (x[row] + x[row + 1])
                          + (x[row + 2] + x[row + 3]);
            sum += row_sum;
        }
    }

    output[index] = (sum * (1.0f / 64.0f) - mean[c]) * scale[c] + bias[c];
}

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() == at::kFloat,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.is_contiguous() && x.dim() == 5,
                "Expected contiguous NCDHW input");

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

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

    TORCH_CHECK(depth >= 4 && height >= 4 && width >= 4,
                "Input is too small for two pooling layers");

    int od = depth / 4;
    int oh = height / 4;
    int ow = width / 4;
    int64_t spatial = int64_t(depth) * height * width;
    int64_t samples = int64_t(batch) * spatial;

    auto output = torch::empty({batch, channels, od, oh, ow}, x.options());
    auto stats = torch::empty({2, channels}, x.options());
    float* mean = stats.data_ptr<float>();
    float* scale = mean + channels;

    if (training) {
        TORCH_CHECK(samples > 1, "Batch normalization needs more than one value");
        int parts = int((samples + CHUNK - 1) / CHUNK);
        auto partial = torch::empty({2, channels, parts}, x.options());
        float* partial_mean = partial.data_ptr<float>();
        float* partial_m2 = partial_mean + int64_t(channels) * parts;

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

        finalize_stats_kernel<<<channels, THREADS, 0, stream>>>(
            partial_mean, partial_m2, weight.data_ptr<float>(),
            running_mean.numel() ? running_mean.data_ptr<float>() : nullptr,
            running_var.numel() ? running_var.data_ptr<float>() : nullptr,
            mean, scale, samples, parts, float(momentum), float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    } else {
        eval_stats_kernel<<<(channels + THREADS - 1) / THREADS,
                            THREADS, 0, stream>>>(
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            weight.data_ptr<float>(), mean, scale, channels, float(eps));
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }

    int64_t total = output.numel();
    if (total > 0) {
        fused_pool_kernel<<<(total + THREADS - 1) / THREADS,
                            THREADS, 0, stream>>>(
            x.data_ptr<float>(), mean, scale, bias.data_ptr<float>(),
            output.data_ptr<float>(), channels, depth, height, width,
            od, oh, ow, total);
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }
    return output;
}
"""

_extension = load_inline(
    name="fused_bn_double_avgpool3d",
    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 bn.weight.dtype != torch.float32
            or min(x.shape[2:]) < 4
            or x.numel() == 0
        ):
            return self.avg_pool2(self.avg_pool1(bn(x)))

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

        use_batch_stats = bn.training or bn.running_mean is None
        running_mean = bn.running_mean
        running_var = bn.running_var
        if running_mean is None:
            running_mean = x.new_empty(0)
            running_var = x.new_empty(0)

        return _extension.bn_pool_cuda(
            x.contiguous(),
            bn.weight,
            bn.bias,
            running_mean,
            running_var,
            use_batch_stats,
            momentum,
            bn.eps,
        )
```