```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_center_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    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 <algorithm>
#include <cstdint>

constexpr int THREADS = 256;
constexpr int CHUNK = 4096;

__device__ __forceinline__ float warp_sum(float value) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        value += __shfl_down_sync(0xffffffff, value, offset);
    return value;
}

template<bool TRAINING>
__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial,
    int spatial,
    int parts)
{
    const int part = blockIdx.x % parts;
    const int nc = blockIdx.x / parts;
    const int begin = part * CHUNK;
    const int end = min(begin + CHUNK, spatial);
    const int64_t base = static_cast<int64_t>(nc) * spatial;

    float sum = 0.0f;
    float square_sum = 0.0f;

    for (int i = begin + threadIdx.x; i < end; i += THREADS) {
        const float value = x[base + i];
        sum += value;
        if (TRAINING)
            square_sum = fmaf(value, value, square_sum);
    }

    sum = warp_sum(sum);
    if (TRAINING)
        square_sum = warp_sum(square_sum);

    __shared__ float sums[8];
    __shared__ float square_sums[8];

    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    if (lane == 0) {
        sums[warp] = sum;
        if (TRAINING)
            square_sums[warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = warp_sum(lane < 8 ? sums[lane] : 0.0f);
        if (TRAINING)
            square_sum = warp_sum(lane < 8 ? square_sums[lane] : 0.0f);

        if (lane == 0) {
            const int64_t offset = static_cast<int64_t>(blockIdx.x) * 2;
            partial[offset] = sum;
            if (TRAINING)
                partial[offset + 1] = square_sum;
        }
    }
}

template<bool TRAINING>
__global__ void finalize_moments_kernel(
    const float* __restrict__ partial,
    const float* __restrict__ weight,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ instance_mean,
    float* __restrict__ scale,
    int batch,
    int channels,
    int spatial,
    int parts,
    float momentum,
    float eps)
{
    const int channel = blockIdx.x;
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    float total_sum = 0.0f;
    float total_square_sum = 0.0f;

    for (int n = warp; n < batch; n += 8) {
        const int nc = n * channels + channel;
        const int64_t offset = static_cast<int64_t>(nc) * parts * 2;

        float sum = 0.0f;
        float square_sum = 0.0f;

        for (int part = lane; part < parts; part += 32) {
            sum += partial[offset + 2 * part];
            if (TRAINING)
                square_sum += partial[offset + 2 * part + 1];
        }

        sum = warp_sum(sum);
        if (TRAINING)
            square_sum = warp_sum(square_sum);

        if (lane == 0) {
            instance_mean[nc] = sum / static_cast<float>(spatial);
            if (TRAINING) {
                total_sum += sum;
                total_square_sum += square_sum;
            }
        }
    }

    if (TRAINING) {
        __shared__ float sums[8];
        __shared__ float square_sums[8];

        if (lane == 0) {
            sums[warp] = total_sum;
            square_sums[warp] = total_square_sum;
        }
        __syncthreads();

        if (warp == 0) {
            float sum = warp_sum(lane < 8 ? sums[lane] : 0.0f);
            float square_sum =
                warp_sum(lane < 8 ? square_sums[lane] : 0.0f);

            if (lane == 0) {
                const double count =
                    static_cast<double>(batch) * spatial;
                const double mean = static_cast<double>(sum) / count;
                const float variance = static_cast<float>(fmax(
                    static_cast<double>(square_sum) / count - mean * mean,
                    0.0));

                scale[channel] = weight[channel] * rsqrtf(variance + eps);

                const float unbiased_variance = static_cast<float>(
                    static_cast<double>(variance) * count / (count - 1.0));

                running_mean[channel] =
                    (1.0f - momentum) * running_mean[channel]
                    + momentum * static_cast<float>(mean);
                running_var[channel] =
                    (1.0f - momentum) * running_var[channel]
                    + momentum * unbiased_variance;
            }
        }
    } else if (threadIdx.x == 0) {
        scale[channel] =
            weight[channel] * rsqrtf(running_var[channel] + eps);
    }
}

__global__ void center_and_scale_kernel(
    const float* __restrict__ x,
    const float* __restrict__ instance_mean,
    const float* __restrict__ scale,
    float* __restrict__ output,
    int spatial,
    int channels)
{
    const int nc = blockIdx.y;
    const float mean = instance_mean[nc];
    const float multiplier = scale[nc % channels];
    const int64_t base = static_cast<int64_t>(nc) * spatial;
    const int begin = blockIdx.x * THREADS * 4 + threadIdx.x;

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        const int i = begin + j * THREADS;
        if (i < spatial)
            output[base + i] = (x[base + i] - mean) * multiplier;
    }
}

torch::Tensor fused_bn_center_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    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 int batch = x.size(0);
    const int channels = x.size(1);
    const int64_t spatial64 = x.size(2) * x.size(3) * x.size(4);

    TORCH_CHECK(batch > 0 && channels > 0 && spatial64 > 0,
                "Expected nonempty input");
    TORCH_CHECK(spatial64 <= INT32_MAX, "Spatial extent is too large");
    TORCH_CHECK(static_cast<int64_t>(batch) * channels <= 65535,
                "Too many batch/channel groups");
    TORCH_CHECK(!training || batch * spatial64 > 1,
                "Expected more than one value per channel when training");

    const int spatial = static_cast<int>(spatial64);
    const int parts = (spatial - 1) / CHUNK + 1;

    auto partial = torch::empty({batch, channels, parts, 2}, x.options());
    auto means = torch::empty({batch, channels}, x.options());
    auto scales = torch::empty({channels}, x.options());
    auto output = torch::empty_like(x);

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        partial_moments_kernel<true>
            <<<batch * channels * parts, THREADS, 0, stream>>>(
                x.data_ptr<float>(), partial.data_ptr<float>(),
                spatial, parts);

        finalize_moments_kernel<true>
            <<<channels, THREADS, 0, stream>>>(
                partial.data_ptr<float>(), weight.data_ptr<float>(),
                running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
                means.data_ptr<float>(), scales.data_ptr<float>(),
                batch, channels, spatial, parts,
                static_cast<float>(momentum), static_cast<float>(eps));
    } else {
        partial_moments_kernel<false>
            <<<batch * channels * parts, THREADS, 0, stream>>>(
                x.data_ptr<float>(), partial.data_ptr<float>(),
                spatial, parts);

        finalize_moments_kernel<false>
            <<<channels, THREADS, 0, stream>>>(
                partial.data_ptr<float>(), weight.data_ptr<float>(),
                running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
                means.data_ptr<float>(), scales.data_ptr<float>(),
                batch, channels, spatial, parts,
                static_cast<float>(momentum), static_cast<float>(eps));
    }

    const dim3 grid(
        (spatial - 1) / (THREADS * 4) + 1,
        batch * channels);

    center_and_scale_kernel<<<grid, THREADS, 0, stream>>>(
        x.data_ptr<float>(), means.data_ptr<float>(),
        scales.data_ptr<float>(), output.data_ptr<float>(),
        spatial, channels);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="conv_transpose3d_fused_bn_center_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_bn_center_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=True
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels,
            out_channels,
            kernel_size,
            stride=stride,
            padding=padding,
            bias=bias,
        )
        self.batch_norm = nn.BatchNorm3d(out_channels)

    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
            or bn.running_mean.dtype != torch.float32
            or x.numel() == 0
            or x.shape[0] * x.shape[1] > 65535
        ):
            x = bn(x)
            return x - x.mean(dim=(2, 3, 4), keepdim=True)

        x = x.contiguous()
        momentum = bn.momentum if bn.momentum is not None else 0.0

        if bn.training:
            if x.numel() // x.shape[1] <= 1:
                x = bn(x)
                return x - x.mean(dim=(2, 3, 4), keepdim=True)

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

        # Spatial centering cancels BatchNorm's channelwise shift:
        # BN(x) - spatial_mean(BN(x))
        # = gamma * inv_std * (x - spatial_mean(x)).
        return _fused_ops.fused_bn_center_cuda(
            x,
            bn.weight,
            bn.running_mean,
            bn.running_var,
            bn.training,
            momentum,
            bn.eps,
        )
```