```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 relu_group_norm_cuda(
    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.h>
#include <cuda_runtime.h>

constexpr int THREADS = 256;
constexpr int REDUCE_TILE = 8192;
constexpr int OUTPUT_TILE = 2048;

__device__ __forceinline__ void warp_reduce(float& a, float& b) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        a += __shfl_down_sync(0xffffffff, a, offset);
        b += __shfl_down_sync(0xffffffff, b, offset);
    }
}

__device__ __forceinline__ void block_reduce(float& a, float& b) {
    __shared__ float sums[8];
    __shared__ float squares[8];

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

    warp_reduce(a, b);
    if (lane == 0) {
        sums[warp] = a;
        squares[warp] = b;
    }
    __syncthreads();

    if (warp == 0) {
        a = lane < 8 ? sums[lane] : 0.0f;
        b = lane < 8 ? squares[lane] : 0.0f;
        warp_reduce(a, b);
    }
}

template <typename scalar_t>
__global__ void relu_partial_moments(
    const scalar_t* __restrict__ x,
    float* __restrict__ partial,
    int64_t group_size,
    int chunks)
{
    const int chunk = blockIdx.x % chunks;
    const int group = blockIdx.x / chunks;
    const int64_t start = static_cast<int64_t>(chunk) * REDUCE_TILE;
    const int64_t base = static_cast<int64_t>(group) * group_size;

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

    #pragma unroll
    for (int j = 0; j < REDUCE_TILE / THREADS; ++j) {
        const int64_t index = start + threadIdx.x + j * THREADS;
        if (index < group_size) {
            float value = static_cast<float>(x[base + index]);
            value = value < 0.0f ? 0.0f : value;
            sum += value;
            square_sum = fmaf(value, value, square_sum);
        }
    }

    block_reduce(sum, square_sum);
    if (threadIdx.x == 0) {
        const int64_t index =
            (static_cast<int64_t>(group) * chunks + chunk) * 2;
        partial[index] = sum;
        partial[index + 1] = square_sum;
    }
}

__global__ void finish_moments(
    const float* __restrict__ partial,
    float* __restrict__ moments,
    int64_t group_size,
    int chunks,
    float eps)
{
    const int group = blockIdx.x;
    const int64_t base = static_cast<int64_t>(group) * chunks * 2;

    float sum = 0.0f;
    float square_sum = 0.0f;
    for (int i = threadIdx.x; i < chunks; i += THREADS) {
        sum += partial[base + i * 2];
        square_sum += partial[base + i * 2 + 1];
    }

    block_reduce(sum, square_sum);
    if (threadIdx.x == 0) {
        const float mean = sum / static_cast<float>(group_size);
        const float variance = fmaxf(
            square_sum / static_cast<float>(group_size) - mean * mean,
            0.0f);
        moments[group * 2] = mean;
        moments[group * 2 + 1] = rsqrtf(variance + eps);
    }
}

template <typename scalar_t>
__global__ void relu_normalize_affine(
    const scalar_t* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ moments,
    scalar_t* __restrict__ output,
    int64_t spatial,
    int channels,
    int channels_per_group,
    int groups,
    int tiles_per_channel)
{
    const int row = blockIdx.x / tiles_per_channel;
    const int tile = blockIdx.x % tiles_per_channel;
    const int channel = row % channels;
    const int batch = row / channels;
    const int group = batch * groups + channel / channels_per_group;

    const float mean = moments[group * 2];
    const float scale = moments[group * 2 + 1] * weight[channel];
    const float shift = bias[channel];
    const int64_t base = static_cast<int64_t>(row) * spatial;
    const int64_t start = static_cast<int64_t>(tile) * OUTPUT_TILE;

    #pragma unroll
    for (int j = 0; j < OUTPUT_TILE / THREADS; ++j) {
        const int64_t index = start + threadIdx.x + j * THREADS;
        if (index < spatial) {
            float value = static_cast<float>(x[base + index]);
            value = value < 0.0f ? 0.0f : value;
            output[base + index] =
                static_cast<scalar_t>(fmaf(value - mean, scale, shift));
        }
    }
}

torch::Tensor relu_group_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps)
{
    TORCH_CHECK(x.is_cuda() && x.is_contiguous(),
                "x must be a contiguous CUDA tensor");
    TORCH_CHECK(x.dim() == 5, "expected NCDHW input");
    TORCH_CHECK(groups > 0 && x.size(1) % groups == 0,
                "invalid group count");
    TORCH_CHECK(weight.device() == x.device() &&
                bias.device() == x.device(),
                "parameters must be on the input device");

    const c10::cuda::CUDAGuard guard(x.device());
    auto output = torch::empty_like(x);
    if (x.numel() == 0) {
        return output;
    }

    const int channels = static_cast<int>(x.size(1));
    const int channels_per_group = channels / groups;
    const int64_t spatial = x.size(2) * x.size(3) * x.size(4);
    const int64_t group_size = spatial * channels_per_group;
    const int chunks = (group_size + REDUCE_TILE - 1) / REDUCE_TILE;
    const int total_groups = x.size(0) * groups;
    const int tiles = (spatial + OUTPUT_TILE - 1) / OUTPUT_TILE;

    auto options = x.options().dtype(torch::kFloat32);
    auto partial = torch::empty({total_groups, chunks, 2}, options);
    auto moments = torch::empty({total_groups, 2}, options);
    auto w = weight.to(torch::kFloat32).contiguous();
    auto b = bias.to(torch::kFloat32).contiguous();

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

    AT_DISPATCH_FLOATING_TYPES_AND2(
        at::ScalarType::Half,
        at::ScalarType::BFloat16,
        x.scalar_type(),
        "relu_group_norm_cuda",
        [&] {
            relu_partial_moments<scalar_t>
                <<<total_groups * chunks, THREADS, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    partial.data_ptr<float>(),
                    group_size,
                    chunks);
            C10_CUDA_KERNEL_LAUNCH_CHECK();

            finish_moments<<<total_groups, THREADS, 0, stream>>>(
                partial.data_ptr<float>(),
                moments.data_ptr<float>(),
                group_size,
                chunks,
                static_cast<float>(eps));
            C10_CUDA_KERNEL_LAUNCH_CHECK();

            relu_normalize_affine<scalar_t>
                <<<x.size(0) * channels * tiles, THREADS, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    w.data_ptr<float>(),
                    b.data_ptr<float>(),
                    moments.data_ptr<float>(),
                    output.data_ptr<scalar_t>(),
                    spatial,
                    channels,
                    channels_per_group,
                    static_cast<int>(groups),
                    tiles);
            C10_CUDA_KERNEL_LAUNCH_CHECK();
        });

    return output;
}
"""

_fused_ops = load_inline(
    name="fused_relu_group_norm_3d_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["relu_group_norm_cuda"],
    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
        )
        self.relu = nn.ReLU()
        self.group_norm = nn.GroupNorm(
            num_groups=groups, num_channels=out_channels
        )

    def forward(self, x):
        x = self.conv_transpose(x)
        if not x.is_cuda or x.dtype not in (
            torch.float32, torch.float16, torch.bfloat16
        ):
            return self.group_norm(F.relu(x))

        return _fused_ops.relu_group_norm_cuda(
            x.contiguous(),
            self.group_norm.weight,
            self.group_norm.bias,
            self.group_norm.num_groups,
            self.group_norm.eps,
        )
```