```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_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    torch::Tensor num_batches_tracked,
    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 <stdint.h>

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;

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

__device__ __forceinline__ Moments block_reduce(
    Moments a, float* means, float* m2s, int* counts) {
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;
    a = warp_reduce(a);

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

    Moments result = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8) {
            result.mean = means[lane];
            result.m2 = m2s[lane];
            result.count = counts[lane];
        }
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial,
    int channels,
    int64_t spatial,
    int chunks,
    int partials_per_channel) {
    constexpr int TILE = 4096;
    const int part = blockIdx.x;
    const int channel = blockIdx.y;
    const int batch = part / chunks;
    const int chunk = part - batch * chunks;
    const int64_t start = int64_t(chunk) * TILE;

    const float* base =
        x + (int64_t(batch) * channels + channel) * spatial;

    Moments local = {0.0f, 0.0f, 0};
    #pragma unroll
    for (int i = 0; i < TILE / 256; ++i) {
        const int64_t position = start + threadIdx.x + i * 256;
        if (position < spatial) {
            const float value = base[position];
            ++local.count;
            const float delta = value - local.mean;
            local.mean += delta / float(local.count);
            local.m2 += delta * (value - local.mean);
        }
    }

    __shared__ float means[8];
    __shared__ float m2s[8];
    __shared__ int counts[8];
    Moments result = block_reduce(local, means, m2s, counts);

    if (threadIdx.x == 0) {
        const int64_t index =
            (int64_t(channel) * partials_per_channel + part) * 3;
        partial[index] = result.mean;
        partial[index + 1] = result.m2;
        partial[index + 2] = float(result.count);
    }
}

__global__ void finalize_moments_kernel(
    const float* __restrict__ partial,
    float* __restrict__ statistics,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    int64_t* __restrict__ num_batches_tracked,
    int partials_per_channel,
    float momentum,
    float eps) {
    const int channel = blockIdx.x;
    Moments local = {0.0f, 0.0f, 0};

    for (int part = threadIdx.x;
         part < partials_per_channel;
         part += blockDim.x) {
        const int64_t index =
            (int64_t(channel) * partials_per_channel + part) * 3;
        Moments other = {
            partial[index],
            partial[index + 1],
            int(partial[index + 2])
        };
        local = merge_moments(local, other);
    }

    __shared__ float means[8];
    __shared__ float m2s[8];
    __shared__ int counts[8];
    Moments result = block_reduce(local, means, m2s, counts);

    if (threadIdx.x == 0) {
        const float variance = fmaxf(result.m2 / float(result.count), 0.0f);
        statistics[channel * 2] = result.mean;
        statistics[channel * 2 + 1] = rsqrtf(variance + eps);

        // A negative momentum denotes cumulative averaging.
        const float factor = momentum < 0.0f
            ? 1.0f / float(*num_batches_tracked + 1)
            : momentum;

        running_mean[channel] =
            (1.0f - factor) * running_mean[channel]
            + factor * result.mean;

        const float unbiased_variance =
            result.m2 / float(result.count - 1);
        running_var[channel] =
            (1.0f - factor) * running_var[channel]
            + factor * unbiased_variance;
    }
}

__global__ void increment_batches_kernel(int64_t* counter) {
    ++*counter;
}

template<bool Training>
__global__ void normalize_pool_kernel(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ statistics,
    const float* __restrict__ running_mean,
    const float* __restrict__ running_var,
    float* __restrict__ output,
    int channels,
    int depth,
    int height,
    int width,
    int output_depth,
    int output_height,
    int output_width,
    int64_t output_elements,
    float eps) {
    const int lane = threadIdx.x & 31;
    const int64_t index =
        int64_t(blockIdx.x) * (blockDim.x / 32) + threadIdx.x / 32;
    if (index >= output_elements) return;

    int64_t remaining = index;
    const int ow = remaining % output_width;
    remaining /= output_width;
    const int oh = remaining % output_height;
    remaining /= output_height;
    const int od = remaining % output_depth;
    remaining /= output_depth;

    const int channel = remaining % channels;
    const int64_t spatial = int64_t(depth) * height * width;

    // Two consecutive stride-2 average pools equal one stride-4
    // average over a 4x4x4 region, including odd-sized inputs.
    const int dz = lane >> 4;
    const int dy = (lane >> 2) & 3;
    const int dx = lane & 3;
    const int64_t position =
        remaining * spatial
        + int64_t(od * 4 + dz) * height * width
        + int64_t(oh * 4 + dy) * width
        + ow * 4 + dx;

    float sum = x[position] + x[position + int64_t(2) * height * width];

    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        sum += __shfl_down_sync(0xffffffff, sum, offset);

    if (lane == 0) {
        float mean;
        float inverse_std;
        if (Training) {
            mean = statistics[channel * 2];
            inverse_std = statistics[channel * 2 + 1];
        } else {
            mean = running_mean[channel];
            inverse_std = rsqrtf(running_var[channel] + eps);
        }

        output[index] =
            ((sum * (1.0f / 64.0f) - mean) * inverse_std)
            * weight[channel] + bias[channel];
    }
}

torch::Tensor bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    torch::Tensor num_batches_tracked,
    bool training,
    double momentum,
    double eps) {
    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 int 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);
    TORCH_CHECK(depth >= 4 && height >= 4 && width >= 4,
                "Input is too small for two average pooling layers");

    const int od = depth / 4;
    const int oh = height / 4;
    const int ow = width / 4;
    const int64_t spatial = int64_t(depth) * height * width;
    TORCH_CHECK(!training || int64_t(batch) * spatial > 1,
                "Batch normalization needs more than one value per channel");

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

    torch::Tensor partial;
    torch::Tensor statistics;
    const float* statistics_ptr = nullptr;

    if (training) {
        const int chunks = int((spatial + 4095) / 4096);
        const int partials_per_channel = batch * chunks;
        partial = torch::empty(
            {channels, partials_per_channel, 3}, x.options());
        statistics = torch::empty({channels, 2}, x.options());
        statistics_ptr = statistics.data_ptr<float>();

        partial_moments_kernel<<<
            dim3(partials_per_channel, channels), 256, 0, stream>>>(
                x.data_ptr<float>(),
                partial.data_ptr<float>(),
                channels,
                spatial,
                chunks,
                partials_per_channel);

        finalize_moments_kernel<<<channels, 256, 0, stream>>>(
            partial.data_ptr<float>(),
            statistics.data_ptr<float>(),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            num_batches_tracked.data_ptr<int64_t>(),
            partials_per_channel,
            float(momentum),
            float(eps));

        increment_batches_kernel<<<1, 1, 0, stream>>>(
            num_batches_tracked.data_ptr<int64_t>());
    }

    const int64_t elements = output.numel();
    const int blocks = int((elements + 7) / 8);

    if (training) {
        normalize_pool_kernel<true><<<blocks, 256, 0, stream>>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            statistics_ptr,
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, depth, height, width,
            od, oh, ow, elements, float(eps));
    } else {
        normalize_pool_kernel<false><<<blocks, 256, 0, stream>>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            nullptr,
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            output.data_ptr<float>(),
            channels, depth, height, width,
            od, oh, ow, elements, float(eps));
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_bn_pool_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 not x.is_contiguous()
            or not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean.dtype != torch.float32
        ):
            return self.avg_pool2(self.avg_pool1(bn(x)))

        return _bn_pool_extension.bn_pool_cuda(
            x,
            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,
        )
```