```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 scale,
    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>

constexpr int THREADS = 256;
constexpr int TILE = 2048;

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

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

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

__device__ __forceinline__ Moments block_reduce(
    Moments value, Moments* shared) {
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

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

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

template<bool TRAINING>
__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int spatial,
    int tiles,
    float scale) {
    const int tile = blockIdx.x % tiles;
    const int sample_channel = blockIdx.x / tiles;
    const int start = tile * TILE;
    const int length = min(TILE, spatial - start);

    const float* values =
        x + static_cast<int64_t>(sample_channel) * spatial + start;

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

    if constexpr (TRAINING) {
        for (int i = threadIdx.x; i < length; i += THREADS) {
            const float value = values[i] * scale;
            ++local.count;
            const float delta = value - local.mean;
            local.mean += delta / float(local.count);
            local.m2 += delta * (value - local.mean);
        }
    } else {
        float sum = 0.0f;
        for (int i = threadIdx.x; i < length; i += THREADS) {
            sum += values[i] * scale;
            ++local.count;
        }
        if (local.count > 0) local.mean = sum / float(local.count);
    }

    __shared__ Moments shared[THREADS / 32];
    const Moments total = block_reduce(local, shared);

    if (threadIdx.x == 0) {
        partial[blockIdx.x] = make_float2(total.mean, total.m2);
    }
}

template<bool TRAINING>
__global__ void finish_bn_pool_kernel(
    const float2* __restrict__ partial,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ output,
    int batch,
    int channels,
    int spatial,
    int tiles,
    float eps,
    float momentum,
    bool update_running) {
    const int channel = blockIdx.x;
    Moments local = {0.0f, 0.0f, 0};

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        Moments sample = {0.0f, 0.0f, 0};
        const int64_t base =
            (static_cast<int64_t>(n) * channels + channel) * tiles;

        for (int tile = 0; tile < tiles; ++tile) {
            const float2 value = partial[base + tile];
            Moments part = {
                value.x,
                value.y,
                min(TILE, spatial - tile * TILE)
            };
            sample = merge_moments(sample, part);
        }

        output[static_cast<int64_t>(n) * channels + channel] = sample.mean;
        if constexpr (TRAINING) {
            local = merge_moments(local, sample);
        }
    }

    __shared__ Moments shared[THREADS / 32];
    __shared__ float normalization_mean;
    __shared__ float normalization_gain;
    __shared__ float normalization_bias;

    if constexpr (TRAINING) {
        const Moments total = block_reduce(local, shared);

        if (threadIdx.x == 0) {
            const float variance =
                fmaxf(total.m2 / float(total.count), 0.0f);

            normalization_mean = total.mean;
            normalization_gain = weight[channel] * rsqrtf(variance + eps);
            normalization_bias = bias[channel];

            if (update_running) {
                const float unbiased_variance =
                    fmaxf(total.m2 / float(total.count - 1), 0.0f);
                running_mean[channel] =
                    (1.0f - momentum) * running_mean[channel]
                    + momentum * total.mean;
                running_var[channel] =
                    (1.0f - momentum) * running_var[channel]
                    + momentum * unbiased_variance;
            }
        }
    } else {
        if (threadIdx.x == 0) {
            normalization_mean = running_mean[channel];
            normalization_gain =
                weight[channel] * rsqrtf(running_var[channel] + eps);
            normalization_bias = bias[channel];
        }
    }

    __syncthreads();

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        const int64_t index =
            static_cast<int64_t>(n) * channels + channel;
        output[index] =
            (output[index] - normalization_mean) * normalization_gain
            + normalization_bias;
    }
}

torch::Tensor fused_bn_pool_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double scale,
    double eps,
    double momentum,
    bool training) {
    const c10::cuda::CUDAGuard guard(x.device());

    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat);
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous());

    const int batch = x.size(0);
    const int channels = x.size(1);
    const int spatial = x.size(2) * x.size(3) * x.size(4);
    const int tiles = (spatial + TILE - 1) / TILE;

    auto partial = torch::empty({batch, channels, tiles, 2}, x.options());
    auto output = torch::empty({batch, channels, 1, 1, 1}, x.options());

    const bool update_running =
        running_mean.numel() != 0 && running_var.numel() != 0;

    float* mean_ptr =
        update_running ? running_mean.data_ptr<float>() : nullptr;
    float* var_ptr =
        update_running ? running_var.data_ptr<float>() : nullptr;

    const auto stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        partial_moments_kernel<true>
            <<<batch * channels * tiles, THREADS, 0, stream>>>(
                x.data_ptr<float>(),
                reinterpret_cast<float2*>(partial.data_ptr<float>()),
                spatial, tiles, static_cast<float>(scale));

        finish_bn_pool_kernel<true>
            <<<channels, THREADS, 0, stream>>>(
                reinterpret_cast<const float2*>(partial.data_ptr<float>()),
                weight.data_ptr<float>(), bias.data_ptr<float>(),
                mean_ptr, var_ptr, output.data_ptr<float>(),
                batch, channels, spatial, tiles,
                static_cast<float>(eps), static_cast<float>(momentum),
                update_running);
    } else {
        TORCH_CHECK(update_running, "Evaluation requires running statistics.");

        partial_moments_kernel<false>
            <<<batch * channels * tiles, THREADS, 0, stream>>>(
                x.data_ptr<float>(),
                reinterpret_cast<float2*>(partial.data_ptr<float>()),
                spatial, tiles, static_cast<float>(scale));

        finish_bn_pool_kernel<false>
            <<<channels, THREADS, 0, stream>>>(
                reinterpret_cast<const float2*>(partial.data_ptr<float>()),
                weight.data_ptr<float>(), bias.data_ptr<float>(),
                mean_ptr, var_ptr, output.data_ptr<float>(),
                batch, channels, spatial, tiles,
                static_cast<float>(eps), static_cast<float>(momentum),
                false);
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

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


class ModelNew(nn.Module):
    def __init__(
        self,
        in_channels,
        out_channels,
        kernel_size,
        scale_factor,
        eps=1e-5,
        momentum=0.1,
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size
        )
        self.scale_factor = scale_factor
        self.batch_norm = nn.BatchNorm3d(
            out_channels, eps=eps, momentum=momentum
        )
        self.global_avg_pool = nn.AdaptiveAvgPool3d((1, 1, 1))

    def forward(self, x):
        x = self.conv_transpose(x)
        bn = self.batch_norm

        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or x.numel() == 0
            or x.numel() // x.shape[1] <= 1
            or not bn.affine
            or bn.weight.dtype != torch.float32
        ):
            return self.global_avg_pool(bn(x * self.scale_factor))

        use_batch_stats = bn.training or (
            bn.running_mean is None and bn.running_var is None
        )

        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())

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

        return _fused_ops.fused_bn_pool_cuda(
            x.contiguous(),
            bn.weight,
            bn.bias,
            running_mean,
            running_var,
            float(self.scale_factor),
            float(bn.eps),
            float(momentum),
            use_batch_stats,
        )
```