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


_cpp_source = """
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_runtime.h>

constexpr int THREADS = 256;
constexpr int STATS_TILE = 4096;
constexpr int OUTPUT_TILE = 1024;

__device__ __forceinline__ void block_reduce_pair(
    float& sum, float& square_sum,
    float* shared_sum, float* shared_square_sum) {
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;

    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        sum += __shfl_down_sync(0xffffffff, sum, offset);
        square_sum += __shfl_down_sync(0xffffffff, square_sum, offset);
    }

    if (lane == 0) {
        shared_sum[warp] = sum;
        shared_square_sum[warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = lane < 8 ? shared_sum[lane] : 0.0f;
        square_sum = lane < 8 ? shared_square_sum[lane] : 0.0f;

        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1) {
            sum += __shfl_down_sync(0xffffffff, sum, offset);
            square_sum += __shfl_down_sync(
                0xffffffff, square_sum, offset);
        }
    }
}

template <typename scalar_t>
__global__ void relu_partial_stats(
    const scalar_t* __restrict__ x,
    float2* __restrict__ partial,
    int64_t group_size,
    int tiles) {
    __shared__ float shared_sum[8];
    __shared__ float shared_square_sum[8];

    const int group = blockIdx.y;
    const int tile = blockIdx.x;
    const int64_t base = int64_t(group) * group_size;
    const int64_t start = int64_t(tile) * STATS_TILE;

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

    #pragma unroll
    for (int offset = 0; offset < STATS_TILE; offset += THREADS) {
        const int64_t index = start + threadIdx.x + offset;
        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_pair(sum, square_sum, shared_sum, shared_square_sum);

    if (threadIdx.x == 0) {
        partial[int64_t(group) * tiles + tile] =
            make_float2(sum, square_sum);
    }
}

__global__ void finish_stats(
    const float2* __restrict__ partial,
    float2* __restrict__ stats,
    int64_t group_size,
    int tiles,
    float eps) {
    __shared__ float shared_sum[8];
    __shared__ float shared_square_sum[8];

    const int group = blockIdx.x;
    float sum = 0.0f;
    float square_sum = 0.0f;

    for (int tile = threadIdx.x; tile < tiles; tile += THREADS) {
        const float2 value = partial[int64_t(group) * tiles + tile];
        sum += value.x;
        square_sum += value.y;
    }

    block_reduce_pair(sum, square_sum, shared_sum, shared_square_sum);

    if (threadIdx.x == 0) {
        const float inv_count = 1.0f / float(group_size);
        const float mean = sum * inv_count;
        float variance = square_sum * inv_count - mean * mean;
        variance = variance < 0.0f ? 0.0f : variance;
        stats[group] = make_float2(mean, 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 float2* __restrict__ stats,
    scalar_t* __restrict__ output,
    int64_t group_size,
    int64_t spatial_size,
    int channels_per_group,
    int groups) {
    const int group = blockIdx.y;
    const int channel_base = (group % groups) * channels_per_group;
    const int64_t base = int64_t(group) * group_size;
    const int64_t start = int64_t(blockIdx.x) * OUTPUT_TILE;
    const float2 moments = stats[group];

    #pragma unroll
    for (int offset = 0; offset < OUTPUT_TILE; offset += THREADS) {
        const int64_t index = start + threadIdx.x + offset;
        if (index < group_size) {
            const int channel = channel_base + index / spatial_size;
            float value = static_cast<float>(x[base + index]);
            value = value < 0.0f ? 0.0f : value;
            const float normalized = (value - moments.x) * moments.y;
            output[base + index] = static_cast<scalar_t>(
                normalized * weight[channel] + bias[channel]);
        }
    }
}

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 an NCDHW tensor");
    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");
    TORCH_CHECK(weight.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat &&
                weight.is_contiguous() && bias.is_contiguous(),
                "Parameters must be contiguous float32 tensors");
    TORCH_CHECK(weight.numel() == x.size(1) &&
                bias.numel() == x.size(1),
                "Incorrect affine parameter sizes");

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

    const int64_t spatial_size = x.size(2) * x.size(3) * x.size(4);
    const int channels_per_group = x.size(1) / groups;
    const int64_t group_size = channels_per_group * spatial_size;
    const int64_t group_count = x.size(0) * groups;
    TORCH_CHECK(group_count <= 65535, "Too many batch groups");

    const int tiles = (group_size + STATS_TILE - 1) / STATS_TILE;
    const int output_tiles =
        (group_size + OUTPUT_TILE - 1) / OUTPUT_TILE;

    auto options = x.options().dtype(torch::kFloat32);
    auto partial = torch::empty({group_count, tiles, 2}, options);
    auto stats = torch::empty({group_count, 2}, options);

    auto partial_ptr =
        reinterpret_cast<float2*>(partial.data_ptr<float>());
    auto stats_ptr =
        reinterpret_cast<float2*>(stats.data_ptr<float>());
    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_stats<scalar_t>
                <<<dim3(tiles, group_count), THREADS, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    partial_ptr,
                    group_size,
                    tiles);
            C10_CUDA_KERNEL_LAUNCH_CHECK();

            finish_stats<<<group_count, THREADS, 0, stream>>>(
                partial_ptr,
                stats_ptr,
                group_size,
                tiles,
                static_cast<float>(eps));
            C10_CUDA_KERNEL_LAUNCH_CHECK();

            relu_normalize_affine<scalar_t>
                <<<dim3(output_tiles, group_count),
                   THREADS, 0, stream>>>(
                    x.data_ptr<scalar_t>(),
                    weight.data_ptr<float>(),
                    bias.data_ptr<float>(),
                    stats_ptr,
                    output.data_ptr<scalar_t>(),
                    group_size,
                    spatial_size,
                    channels_per_group,
                    groups);
            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)
            or x.shape[0] * self.group_norm.num_groups > 65535
            or (
                torch.is_grad_enabled()
                and (
                    x.requires_grad
                    or self.group_norm.weight.requires_grad
                    or self.group_norm.bias.requires_grad
                )
            )
        ):
            return self.group_norm(self.relu(x))

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