```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_relu_group_norm(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    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>

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

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

__device__ __forceinline__ void reduce_pair(float& a, float& b) {
    __shared__ float sa[8];
    __shared__ float sb[8];

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

    a = warp_sum(a);
    b = warp_sum(b);
    if (lane == 0) {
        sa[warp] = a;
        sb[warp] = b;
    }
    __syncthreads();

    if (warp == 0) {
        a = lane < 8 ? sa[lane] : 0.0f;
        b = lane < 8 ? sb[lane] : 0.0f;
        a = warp_sum(a);
        b = warp_sum(b);
    }
}

__device__ __forceinline__ float4 relu4(float4 v) {
    v.x = v.x < 0.0f ? 0.0f : v.x;
    v.y = v.y < 0.0f ? 0.0f : v.y;
    v.z = v.z < 0.0f ? 0.0f : v.z;
    v.w = v.w < 0.0f ? 0.0f : v.w;
    return v;
}

// Input is physically NDHWC. Each block reduces a spatial chunk
// of one normalization group.
__global__ void partial_stats_kernel(
    const float4* __restrict__ x,
    float2* __restrict__ partial,
    int spatial,
    int channel_vectors,
    int group_vectors,
    int groups,
    int chunks)
{
    const int ng = blockIdx.y;
    const int n = ng / groups;
    const int g = ng % groups;
    const int chunk = blockIdx.x;
    const int group_length = spatial * group_vectors;

    const int64_t base =
        static_cast<int64_t>(n) * spatial * channel_vectors
        + g * group_vectors;

    float sum = 0.0f;
    float squares = 0.0f;

    #pragma unroll
    for (int j = 0; j < CHUNK / THREADS; ++j) {
        const int i = chunk * CHUNK + threadIdx.x + j * THREADS;
        if (i < group_length) {
            const int s = i / group_vectors;
            const int c = i % group_vectors;
            const float4 v = relu4(
                x[base + static_cast<int64_t>(s) * channel_vectors + c]);

            sum += (v.x + v.y) + (v.z + v.w);
            squares += (v.x * v.x + v.y * v.y)
                     + (v.z * v.z + v.w * v.w);
        }
    }

    reduce_pair(sum, squares);
    if (threadIdx.x == 0)
        partial[static_cast<int64_t>(ng) * chunks + chunk] =
            make_float2(sum, squares);
}

__global__ void finish_stats_kernel(
    const float2* __restrict__ partial,
    float2* __restrict__ stats,
    int chunks,
    int64_t elements_per_group,
    float eps)
{
    const int ng = blockIdx.x;
    float sum = 0.0f;
    float squares = 0.0f;

    for (int i = threadIdx.x; i < chunks; i += THREADS) {
        const float2 p = partial[static_cast<int64_t>(ng) * chunks + i];
        sum += p.x;
        squares += p.y;
    }

    reduce_pair(sum, squares);
    if (threadIdx.x == 0) {
        const float reciprocal = 1.0f / static_cast<float>(elements_per_group);
        const float mean = sum * reciprocal;
        const float variance = fmaxf(squares * reciprocal - mean * mean, 0.0f);
        stats[ng] = make_float2(mean, rsqrtf(variance + eps));
    }
}

__global__ void normalize_kernel(
    const float4* __restrict__ x,
    float4* __restrict__ output,
    const float2* __restrict__ stats,
    const float4* __restrict__ weight,
    const float4* __restrict__ bias,
    int vectors_per_batch,
    int channel_vectors,
    int group_vectors,
    int groups)
{
    const int n = blockIdx.y;
    const int64_t base = static_cast<int64_t>(n) * vectors_per_batch;

    #pragma unroll
    for (int j = 0; j < CHUNK / THREADS; ++j) {
        const int i = blockIdx.x * CHUNK + threadIdx.x + j * THREADS;
        if (i < vectors_per_batch) {
            const int c = i % channel_vectors;
            const int g = c / group_vectors;

            const float2 st = stats[n * groups + g];
            const float4 w = weight[c];
            const float4 b = bias[c];
            float4 v = relu4(x[base + i]);

            v.x = (v.x - st.x) * (st.y * w.x) + b.x;
            v.y = (v.y - st.x) * (st.y * w.y) + b.y;
            v.z = (v.z - st.x) * (st.y * w.z) + b.z;
            v.w = (v.w - st.x) * (st.y * w.w) + b.w;

            output[base + i] = v;
        }
    }
}

torch::Tensor fused_relu_group_norm(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps)
{
    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 int spatial = static_cast<int>(x.size(2) * x.size(3) * x.size(4));
    const int channel_vectors = channels / 4;
    const int group_vectors = channel_vectors / static_cast<int>(groups);
    const int chunks = (spatial * group_vectors + CHUNK - 1) / CHUNK;
    const int total_groups = batch * static_cast<int>(groups);
    const int vectors_per_batch = spatial * channel_vectors;

    auto output = torch::empty_like(x);
    auto partial = torch::empty({total_groups, chunks, 2}, x.options());
    auto stats = torch::empty({total_groups, 2}, x.options());

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

    partial_stats_kernel<<<dim3(chunks, total_groups), THREADS, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        spatial,
        channel_vectors,
        group_vectors,
        static_cast<int>(groups),
        chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_stats_kernel<<<total_groups, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        reinterpret_cast<float2*>(stats.data_ptr<float>()),
        chunks,
        static_cast<int64_t>(spatial) * (channels / groups),
        static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    normalize_kernel<<<
        dim3((vectors_per_batch + CHUNK - 1) / CHUNK, batch),
        THREADS, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<float4*>(output.data_ptr<float>()),
        reinterpret_cast<const float2*>(stats.data_ptr<float>()),
        reinterpret_cast<const float4*>(weight.data_ptr<float>()),
        reinterpret_cast<const float4*>(bias.data_ptr<float>()),
        vectors_per_batch,
        channel_vectors,
        group_vectors,
        static_cast<int>(groups));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_fused_ops = load_inline(
    name="relu_groupnorm_ndhwc_v2",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_relu_group_norm"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_channels, out_channels, kernel_size, groups, bias=False
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size, bias=bias
        ).to(memory_format=torch.channels_last_3d)
        self.relu = nn.ReLU()
        self.group_norm = nn.GroupNorm(groups, out_channels)

    def forward(self, x):
        if x.is_cuda:
            x = x.contiguous(memory_format=torch.channels_last_3d)

        x = self.conv_transpose(x)
        norm = self.group_norm

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.numel() > 0
            and (x.shape[1] // norm.num_groups) % 4 == 0
            and x.shape[0] * norm.num_groups <= 65535
            and x[0].numel() < 2**31
            and norm.weight.dtype == torch.float32
            and norm.bias.dtype == torch.float32
        ):
            x = x.contiguous(memory_format=torch.channels_last_3d)
            return _fused_ops.fused_relu_group_norm(
                x,
                norm.weight.contiguous(),
                norm.bias.contiguous(),
                norm.num_groups,
                norm.eps,
            )

        return F.group_norm(
            F.relu(x), norm.num_groups, norm.weight, norm.bias, norm.eps
        )
```