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


_cpp_source = r"""
torch::Tensor scaled_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>
#include <algorithm>
#include <cstdint>

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

template<bool Training>
__global__ void partial_moments(
    const float* __restrict__ x,
    double2* __restrict__ partial,
    int64_t spatial,
    int tiles,
    float scale)
{
    const int tile = blockIdx.x % tiles;
    const int64_t nc = blockIdx.x / tiles;
    const int64_t begin = int64_t(tile) * TILE;
    const int64_t end = min(begin + int64_t(TILE), spatial);

    double sum = 0.0;
    double square_sum = 0.0;
    for (int64_t s = begin + threadIdx.x; s < end; s += THREADS) {
        const float value = x[nc * spatial + s] * scale;
        const double v = static_cast<double>(value);
        sum += v;
        if (Training) {
            square_sum += v * v;
        }
    }

    __shared__ double sums[THREADS];
    __shared__ double squares[THREADS];
    sums[threadIdx.x] = sum;
    if (Training) {
        squares[threadIdx.x] = square_sum;
    }
    __syncthreads();

    for (int offset = THREADS / 2; offset > 0; offset >>= 1) {
        if (threadIdx.x < offset) {
            sums[threadIdx.x] += sums[threadIdx.x + offset];
            if (Training) {
                squares[threadIdx.x] += squares[threadIdx.x + offset];
            }
        }
        __syncthreads();
    }

    if (threadIdx.x == 0) {
        partial[nc * tiles + tile] = make_double2(
            sums[0], Training ? squares[0] : 0.0);
    }
}

template<bool Training>
__global__ void normalize_and_pool(
    const double2* __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,
    int64_t spatial,
    int tiles,
    float eps,
    float momentum)
{
    const int c = blockIdx.x;
    __shared__ float mean;
    __shared__ float inv_std;
    __shared__ double sums[THREADS];
    __shared__ double squares[THREADS];

    if (Training) {
        double sum = 0.0;
        double square_sum = 0.0;

        for (int q = threadIdx.x; q < batch * tiles; q += THREADS) {
            const int n = q / tiles;
            const int t = q % tiles;
            const double2 p =
                partial[(int64_t(n) * channels + c) * tiles + t];
            sum += p.x;
            square_sum += p.y;
        }

        sums[threadIdx.x] = sum;
        squares[threadIdx.x] = square_sum;
        __syncthreads();

        for (int offset = THREADS / 2; offset > 0; offset >>= 1) {
            if (threadIdx.x < offset) {
                sums[threadIdx.x] += sums[threadIdx.x + offset];
                squares[threadIdx.x] += squares[threadIdx.x + offset];
            }
            __syncthreads();
        }

        if (threadIdx.x == 0) {
            const double count = double(batch) * double(spatial);
            const double m = sums[0] / count;
            const double variance =
                fmax(squares[0] / count - m * m, 0.0);

            mean = static_cast<float>(m);
            inv_std = rsqrtf(static_cast<float>(variance) + eps);

            const float unbiased =
                static_cast<float>(variance * count / (count - 1.0));
            running_mean[c] =
                (1.0f - momentum) * running_mean[c] + momentum * mean;
            running_var[c] =
                (1.0f - momentum) * running_var[c] + momentum * unbiased;
        }
    } else {
        if (threadIdx.x == 0) {
            mean = running_mean[c];
            inv_std = rsqrtf(running_var[c] + eps);
        }
    }
    __syncthreads();

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        double sum = 0.0;
        const int64_t base = (int64_t(n) * channels + c) * tiles;
        for (int t = 0; t < tiles; ++t) {
            sum += partial[base + t].x;
        }

        const double centered_mean = sum / double(spatial) - double(mean);
        output[int64_t(n) * channels + c] = static_cast<float>(
            centered_mean * double(inv_std) * double(weight[c])
            + double(bias[c]));
    }
}

torch::Tensor scaled_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());

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

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

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

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    auto partial_ptr = reinterpret_cast<double2*>(partial.data_ptr<double>());
    const int blocks = batch * channels * tiles;

    if (training) {
        partial_moments<true><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(), partial_ptr, spatial, tiles,
            static_cast<float>(scale));
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        normalize_and_pool<true><<<channels, THREADS, 0, stream>>>(
            partial_ptr,
            weight.data_ptr<float>(), bias.data_ptr<float>(),
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            output.data_ptr<float>(), batch, channels, spatial, tiles,
            static_cast<float>(eps), static_cast<float>(momentum));
    } else {
        partial_moments<false><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(), partial_ptr, spatial, tiles,
            static_cast<float>(scale));
        C10_CUDA_KERNEL_LAUNCH_CHECK();

        normalize_and_pool<false><<<channels, THREADS, 0, stream>>>(
            partial_ptr,
            weight.data_ptr<float>(), bias.data_ptr<float>(),
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            output.data_ptr<float>(), batch, channels, spatial, tiles,
            static_cast<float>(eps), static_cast<float>(momentum));
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="scaled_bn_global_pool_3d",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["scaled_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

        # Preserve autograd and support other devices, dtypes, and layouts.
        use_custom = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and x.numel() > 0
            and bn.affine
            and bn.track_running_stats
            and bn.weight.dtype == torch.float32
            and bn.running_mean.dtype == torch.float32
            and not (
                torch.is_grad_enabled()
                and (
                    x.requires_grad
                    or bn.weight.requires_grad
                    or bn.bias.requires_grad
                )
            )
        )
        if not use_custom:
            return self.global_avg_pool(bn(x * self.scale_factor))

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

        return _fused_ops.scaled_bn_pool_cuda(
            x,
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            float(self.scale_factor),
            float(bn.eps),
            float(momentum),
            bn.training,
        )
```