```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,
    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 CHUNK = 2048;

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

__device__ __forceinline__ Stats combine(Stats a, Stats 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);

    Stats 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__ Stats warp_reduce(Stats value) {
    const int lane = threadIdx.x & 31;
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Stats 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 = combine(value, other);
    }
    return value;
}

__global__ void partial_stats_kernel(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int64_t* batches,
    int spatial,
    int parts,
    float scale,
    bool training)
{
    const int part = blockIdx.x % parts;
    const int sample_channel = blockIdx.x / parts;
    const int begin = part * CHUNK;
    const int end = min(begin + CHUNK, spatial);
    const int64_t base = int64_t(sample_channel) * spatial;

    Stats value = {0.0f, 0.0f, 0};
    for (int i = begin + threadIdx.x; i < end; i += THREADS) {
        const float v = x[base + i] * scale;
        ++value.count;
        const float delta = v - value.mean;
        value.mean += delta / float(value.count);
        value.m2 += delta * (v - value.mean);
    }

    value = warp_reduce(value);

    __shared__ Stats warp_stats[THREADS / 32];
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

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

    if (warp == 0) {
        value = lane < THREADS / 32
            ? warp_stats[lane]
            : Stats{0.0f, 0.0f, 0};
        value = warp_reduce(value);
        if (lane == 0) {
            partial[blockIdx.x] = make_float2(value.mean, value.m2);
        }
    }

    if (training && blockIdx.x == 0 && threadIdx.x == 0) {
        batches[0] += 1;
    }
}

__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,
    const int64_t* batches,
    float* __restrict__ output,
    int batch,
    int channels,
    int spatial,
    int parts,
    float eps,
    float momentum,
    bool training)
{
    const int channel = blockIdx.x;
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    extern __shared__ float shared[];
    float* sample_mean = shared;
    float* sample_m2 = shared + batch;
    float* parameters = shared + 2 * batch;

    for (int n = warp; n < batch; n += THREADS / 32) {
        Stats value = {0.0f, 0.0f, 0};
        const int64_t base = (int64_t(n) * channels + channel) * parts;

        for (int p = lane; p < parts; p += 32) {
            const float2 v = partial[base + p];
            const int count = min(CHUNK, spatial - p * CHUNK);
            value = combine(value, Stats{v.x, v.y, count});
        }

        value = warp_reduce(value);
        if (lane == 0) {
            sample_mean[n] = value.mean;
            sample_m2[n] = value.m2;
        }
    }
    __syncthreads();

    if (threadIdx.x == 0) {
        float mean;
        float variance;

        if (training) {
            Stats total = {0.0f, 0.0f, 0};
            for (int n = 0; n < batch; ++n) {
                total = combine(
                    total,
                    Stats{sample_mean[n], sample_m2[n], spatial});
            }

            mean = total.mean;
            variance = fmaxf(total.m2 / float(total.count), 0.0f);

            const float factor = momentum < 0.0f
                ? 1.0f / float(batches[0])
                : momentum;
            const float unbiased =
                fmaxf(total.m2 / float(total.count - 1), 0.0f);

            running_mean[channel] =
                (1.0f - factor) * running_mean[channel] + factor * mean;
            running_var[channel] =
                (1.0f - factor) * running_var[channel] + factor * unbiased;
        } else {
            mean = running_mean[channel];
            variance = running_var[channel];
        }

        parameters[0] = mean;
        parameters[1] = rsqrtf(variance + eps) * weight[channel];
        parameters[2] = bias[channel];
    }
    __syncthreads();

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        output[int64_t(n) * channels + channel] =
            (sample_mean[n] - parameters[0]) * parameters[1]
            + parameters[2];
    }
}

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,
    double scale,
    double eps,
    double momentum,
    bool training)
{
    const c10::cuda::CUDAGuard guard(x.device());

    TORCH_CHECK(x.is_contiguous(), "Expected contiguous input");
    TORCH_CHECK(x.scalar_type() == torch::kFloat32,
                "Expected float32 input");

    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 parts = (spatial + CHUNK - 1) / CHUNK;

    TORCH_CHECK(!training || int64_t(batch) * spatial > 1,
                "Batch normalization requires more than one value per channel");

    auto partial = torch::empty({batch, channels, parts, 2}, x.options());
    auto output = torch::empty({batch, channels, 1, 1, 1}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    partial_stats_kernel<<<batch * channels * parts, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        batches.data_ptr<int64_t>(),
        spatial,
        parts,
        static_cast<float>(scale),
        training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_bn_pool_kernel<<<
        channels, THREADS, (2 * batch + 3) * sizeof(float), stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        batches.data_ptr<int64_t>(),
        output.data_ptr<float>(),
        batch,
        channels,
        spatial,
        parts,
        static_cast<float>(eps),
        static_cast<float>(momentum),
        training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_fused_ops = load_inline(
    name="conv3d_transpose_fused_bn_pool_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,
        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 not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
        ):
            return self.global_avg_pool(bn(x * self.scale_factor))

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