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


_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 <cmath>
#include <cstdint>

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

__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;
}

__global__ void partial_stats_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_sum,
    float* __restrict__ partial_square,
    int64_t spatial,
    int chunks,
    int64_t tasks)
{
    __shared__ float shared_sum[8];
    __shared__ float shared_square[8];

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

    for (int64_t task = blockIdx.x; task < tasks; task += gridDim.x) {
        const int64_t row = task / chunks;
        const int chunk = task % chunks;
        const int64_t start = static_cast<int64_t>(chunk) * TILE;
        const int64_t base = row * spatial;

        float sum = 0.0f;
        float square = 0.0f;

        #pragma unroll
        for (int k = 0; k < TILE / THREADS; ++k) {
            const int64_t position = start + threadIdx.x + k * THREADS;
            if (position < spatial) {
                const float value = x[base + position];
                sum += value;
                square = fmaf(value, value, square);
            }
        }

        sum = warp_sum(sum);
        square = warp_sum(square);

        if (lane == 0) {
            shared_sum[warp] = sum;
            shared_square[warp] = square;
        }
        __syncthreads();

        if (warp == 0) {
            sum = lane < 8 ? shared_sum[lane] : 0.0f;
            square = lane < 8 ? shared_square[lane] : 0.0f;
            sum = warp_sum(sum);
            square = warp_sum(square);
            if (lane == 0) {
                partial_sum[task] = sum;
                partial_square[task] = square;
            }
        }
        __syncthreads();
    }
}

__global__ void finalize_stats_kernel(
    const float* __restrict__ partial_sum,
    const float* __restrict__ partial_square,
    const float* __restrict__ weight,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    int64_t* __restrict__ batches,
    float* __restrict__ instance_mean,
    float* __restrict__ scale,
    int batch,
    int channels,
    int64_t spatial,
    int chunks,
    float eps,
    double momentum,
    bool training)
{
    __shared__ float shared_sum[8];
    __shared__ float shared_square[8];

    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 = 0.0f;

    for (int n = warp; n < batch; n += 8) {
        const int64_t row = static_cast<int64_t>(n) * channels + channel;
        const int64_t base = row * chunks;
        float sum = 0.0f;
        float square = 0.0f;

        for (int p = lane; p < chunks; p += 32) {
            sum += partial_sum[base + p];
            square += partial_square[base + p];
        }

        sum = warp_sum(sum);
        square = warp_sum(square);

        if (lane == 0) {
            instance_mean[row] = sum / static_cast<float>(spatial);
            total_sum += sum;
            total_square += square;
        }
    }

    if (lane == 0) {
        shared_sum[warp] = total_sum;
        shared_square[warp] = total_square;
    }
    __syncthreads();

    if (warp == 0) {
        float sum = lane < 8 ? shared_sum[lane] : 0.0f;
        float square = lane < 8 ? shared_square[lane] : 0.0f;
        sum = warp_sum(sum);
        square = warp_sum(square);

        if (lane == 0) {
            float variance = running_var[channel];

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

                const double factor = momentum < 0.0
                    ? 1.0 / static_cast<double>(*batches + 1)
                    : momentum;

                running_mean[channel] = static_cast<float>(
                    (1.0 - factor) * running_mean[channel] + factor * mean);
                running_var[channel] = static_cast<float>(
                    (1.0 - factor) * running_var[channel]
                    + factor * biased_var * count / (count - 1.0));
            }

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

__global__ void normalize_center_kernel(
    const float* __restrict__ x,
    const float* __restrict__ instance_mean,
    const float* __restrict__ scale,
    float* __restrict__ out,
    int64_t* __restrict__ batches,
    int64_t elements,
    int64_t spatial,
    int channels,
    bool training)
{
    const int64_t thread =
        static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
    const int64_t step = static_cast<int64_t>(gridDim.x) * blockDim.x;
    const int64_t vectors = elements / 4;

    if (training && thread == 0)
        *batches += 1;

    for (int64_t v = thread; v < vectors; v += step) {
        const int64_t i = v * 4;
        const float4 values = reinterpret_cast<const float4*>(x)[v];

        const int64_t row0 = i / spatial;
        const int64_t row3 = (i + 3) / spatial;
        float4 result;

        if (row0 == row3) {
            const float mean = instance_mean[row0];
            const float factor = scale[row0 % channels];
            result.x = (values.x - mean) * factor;
            result.y = (values.y - mean) * factor;
            result.z = (values.z - mean) * factor;
            result.w = (values.w - mean) * factor;
        } else {
            const int64_t row1 = (i + 1) / spatial;
            const int64_t row2 = (i + 2) / spatial;
            result.x = (values.x - instance_mean[row0]) * scale[row0 % channels];
            result.y = (values.y - instance_mean[row1]) * scale[row1 % channels];
            result.z = (values.z - instance_mean[row2]) * scale[row2 % channels];
            result.w = (values.w - instance_mean[row3]) * scale[row3 % channels];
        }

        reinterpret_cast<float4*>(out)[v] = result;
    }

    for (int64_t i = vectors * 4 + thread; i < elements; i += step) {
        const int64_t row = i / spatial;
        out[i] = (x[i] - instance_mean[row]) * scale[row % channels];
    }
}

torch::Tensor fused_bn_center_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    torch::Tensor batches,
    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 chunks = static_cast<int>((spatial + TILE - 1) / TILE);
    const int64_t rows = static_cast<int64_t>(batch) * channels;
    const int64_t tasks = rows * chunks;

    auto partial_sum = torch::empty({tasks}, x.options());
    auto partial_square = torch::empty({tasks}, x.options());
    auto instance_mean = torch::empty({rows}, x.options());
    auto scale = torch::empty({channels}, x.options());
    auto out = torch::empty_like(x);

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

    partial_stats_kernel<<<
        static_cast<int>(std::min<int64_t>(tasks, 65535)),
        THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        partial_sum.data_ptr<float>(),
        partial_square.data_ptr<float>(),
        spatial, chunks, tasks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finalize_stats_kernel<<<channels, THREADS, 0, stream>>>(
        partial_sum.data_ptr<float>(),
        partial_square.data_ptr<float>(),
        weight.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        batches.data_ptr<int64_t>(),
        instance_mean.data_ptr<float>(),
        scale.data_ptr<float>(),
        batch, channels, spatial, chunks,
        static_cast<float>(eps), momentum, training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const int blocks = static_cast<int>(std::min<int64_t>(
        (x.numel() + THREADS * 4 - 1) / (THREADS * 4), 65535));

    normalize_center_kernel<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        instance_mean.data_ptr<float>(),
        scale.data_ptr<float>(),
        out.data_ptr<float>(),
        batches.data_ptr<int64_t>(),
        x.numel(), spatial, channels, training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return out;
}
"""

_cpp_source = r"""
torch::Tensor fused_bn_center_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    torch::Tensor batches,
    double eps,
    double momentum,
    bool training);
"""

_fused_ops = load_inline(
    name="fused_bn_spatial_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

        # Preserve autograd behavior when gradients are requested.
        needs_grad = torch.is_grad_enabled() and (
            x.requires_grad
            or (bn.weight is not None and bn.weight.requires_grad)
            or (bn.bias is not None and bn.bias.requires_grad)
        )

        supported = (
            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 is not None
            and bn.weight.dtype == torch.float32
            and bn.weight.is_contiguous()
            and bn.running_mean is not None
            and bn.running_var is not None
            and bn.num_batches_tracked is not None
            and bn.running_mean.dtype == torch.float32
            and bn.running_var.dtype == torch.float32
            and (not bn.training or x.numel() // x.size(1) > 1)
        )

        if needs_grad or not supported:
            x = bn(x)
            return x - x.mean(dim=(2, 3, 4), keepdim=True)

        # Spatial centering cancels the BN mean and affine bias.
        return _fused_ops.fused_bn_center_cuda(
            x,
            bn.weight,
            bn.running_mean,
            bn.running_var,
            bn.num_batches_tracked,
            bn.eps,
            -1.0 if bn.momentum is None else bn.momentum,
            bn.training,
        )
```