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


_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor fused_bias_scale_bn_pool(
    torch::Tensor x,
    torch::Tensor conv_bias,
    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 <cstdint>

constexpr int THREADS = 256;
constexpr int CHUNK = 1024;

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;
    int count = a.count + b.count;
    float delta = b.mean - a.mean;
    float ratio = float(b.count) / float(count);
    return {
        a.mean + delta * ratio,
        a.m2 + b.m2 + delta * delta * float(a.count) * ratio,
        count
    };
}

__device__ __forceinline__ Moments warp_reduce(Moments a) {
    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) {
    __shared__ float means[8];
    __shared__ float m2s[8];
    __shared__ int counts[8];

    int lane = threadIdx.x & 31;
    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();

    a = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8)
            a = {means[lane], m2s[lane], counts[lane]};
        a = warp_reduce(a);
    }
    return a;
}

// Input is physically NDHWC. Adjacent lanes read adjacent channels.
__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    const float* __restrict__ conv_bias,
    float2* __restrict__ partial,
    int channels,
    int spatial,
    int chunks,
    float scale)
{
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;
    const int channel = blockIdx.x * 32 + lane;
    const int chunk = blockIdx.y;
    const int n = blockIdx.z;
    const int begin = chunk * CHUNK;
    const int end = min(begin + CHUNK, spatial);

    __shared__ double sums[8][32];
    __shared__ double squares[8][32];

    double sum = 0.0;
    double square = 0.0;

    if (channel < channels) {
        const float b = conv_bias[channel];
        const int64_t base = int64_t(n) * spatial * channels + channel;

        for (int s = begin + warp; s < end; s += 8) {
            float value = __fmul_rn(
                __fadd_rn(x[base + int64_t(s) * channels], b),
                scale);
            double v = double(value);
            sum += v;
            square = fma(v, v, square);
        }
    }

    sums[warp][lane] = sum;
    squares[warp][lane] = square;
    __syncthreads();

    if (warp == 0 && channel < channels) {
        double total = 0.0;
        double total_square = 0.0;

        #pragma unroll
        for (int w = 0; w < 8; ++w) {
            total += sums[w][lane];
            total_square += squares[w][lane];
        }

        double mean = total / double(end - begin);
        double m2 = fmax(total_square - total * mean, 0.0);

        partial[(int64_t(n) * channels + channel) * chunks + chunk] =
            make_float2(float(mean), float(m2));
    }
}

__global__ void normalize_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 chunks,
    float eps,
    float momentum,
    bool training)
{
    int channel = blockIdx.x;
    __shared__ float mean;
    __shared__ float invstd;

    if (training) {
        Moments total = {0.0f, 0.0f, 0};

        for (int i = threadIdx.x; i < batch * chunks; i += THREADS) {
            int n = i / chunks;
            int chunk = i % chunks;
            int count = min(CHUNK, spatial - chunk * CHUNK);
            float2 p = partial[
                (int64_t(n) * channels + channel) * chunks + chunk];
            total = merge_moments(total, {p.x, p.y, count});
        }

        total = block_reduce(total);

        if (threadIdx.x == 0) {
            mean = total.mean;
            invstd = rsqrtf(fmaxf(total.m2 / float(total.count), 0.0f) + eps);
            running_mean[channel] =
                (1.0f - momentum) * running_mean[channel] +
                momentum * total.mean;
            running_var[channel] =
                (1.0f - momentum) * running_var[channel] +
                momentum * total.m2 / float(total.count - 1);
        }
    } else if (threadIdx.x == 0) {
        mean = running_mean[channel];
        invstd = rsqrtf(running_var[channel] + eps);
    }

    __syncthreads();

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        double sum = 0.0;
        int64_t base = (int64_t(n) * channels + channel) * chunks;

        for (int chunk = 0; chunk < chunks; ++chunk) {
            int count = min(CHUNK, spatial - chunk * CHUNK);
            sum += double(partial[base + chunk].x) * double(count);
        }

        float sample_mean = float(sum / double(spatial));
        output[int64_t(n) * channels + channel] =
            ((sample_mean - mean) * invstd) * weight[channel] +
            bias[channel];
    }
}

torch::Tensor fused_bias_scale_bn_pool(
    torch::Tensor x,
    torch::Tensor conv_bias,
    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 device_guard(x.device());

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

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

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

    dim3 grid((channels + 31) / 32, chunks, batch);
    partial_moments_kernel<<<grid, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        conv_bias.data_ptr<float>(),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        channels,
        spatial,
        chunks,
        float(scale));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    normalize_pool_kernel<<<channels, THREADS, 0, 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>(),
        output.data_ptr<float>(),
        batch,
        channels,
        spatial,
        chunks,
        float(eps),
        float(momentum),
        training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_extension = None


def _get_extension():
    global _extension
    if _extension is None:
        _extension = load_inline(
            name="transpose3d_channels_last_fused_bn_pool_v2",
            cpp_sources=_cpp_source,
            cuda_sources=_cuda_source,
            functions=["fused_bias_scale_bn_pool"],
            extra_cflags=["-O3"],
            extra_cuda_cflags=["-O3"],
            verbose=False,
        )
    return _extension


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.conv_transpose.to(memory_format=torch.channels_last_3d)
        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):
        conv = self.conv_transpose
        bn = self.batch_norm

        tensors = (
            conv.bias,
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
        )
        supported = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 5
            and x.numel() > 0
            and not torch.is_autocast_enabled()
            and conv.weight.dtype == torch.float32
            and all(
                t is not None
                and t.device == x.device
                and t.dtype == torch.float32
                and t.is_contiguous()
                for t in tensors
            )
        )

        if not supported:
            return self.global_avg_pool(
                bn(conv(x) * self.scale_factor)
            )

        x = x.contiguous(memory_format=torch.channels_last_3d)
        weight = conv.weight.contiguous(memory_format=torch.channels_last_3d)

        with torch.backends.cudnn.flags(
            enabled=torch.backends.cudnn.enabled,
            benchmark=True,
            deterministic=torch.backends.cudnn.deterministic,
            allow_tf32=torch.backends.cudnn.allow_tf32,
        ):
            x = F.conv_transpose3d(
                x,
                weight,
                None,
                conv.stride,
                conv.padding,
                conv.output_padding,
                conv.groups,
                conv.dilation,
            )

        x = x.contiguous(memory_format=torch.channels_last_3d)

        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 _get_extension().fused_bias_scale_bn_pool(
            x,
            conv.bias,
            bn.weight,
            bn.bias,
            bn.running_mean,
            bn.running_var,
            float(self.scale_factor),
            bn.eps,
            momentum,
            bn.training,
        )
```