```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_runtime.h>

constexpr int THREADS = 256;
constexpr int STATS_VECTORS = 2048;
constexpr int OUTPUT_VECTORS = 1024;

__device__ __forceinline__ float relu_value(float x) {
    return x < 0.0f ? 0.0f : x;
}

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

__device__ __forceinline__ void block_sum_pair(
    float a, float b, float& total_a, float& total_b) {
    __shared__ float shared_a[8];
    __shared__ float shared_b[8];

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

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

    if (warp == 0) {
        a = lane < 8 ? shared_a[lane] : 0.0f;
        b = lane < 8 ? shared_b[lane] : 0.0f;
        a = warp_sum(a);
        b = warp_sum(b);
    }
    total_a = a;
    total_b = b;
}

__global__ void partial_statistics(
    const float4* __restrict__ x,
    float2* __restrict__ partial,
    int64_t vectors_per_group,
    int chunks) {
    const int group = blockIdx.y;
    const int chunk = blockIdx.x;
    const int64_t start = int64_t(chunk) * STATS_VECTORS;

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

    #pragma unroll
    for (int i = 0; i < STATS_VECTORS / THREADS; ++i) {
        const int64_t index = start + threadIdx.x + i * THREADS;
        if (index < vectors_per_group) {
            const float4 v = x[int64_t(group) * vectors_per_group + index];
            const float a = relu_value(v.x);
            const float b = relu_value(v.y);
            const float c = relu_value(v.z);
            const float d = relu_value(v.w);

            sum += (a + b) + (c + d);
            square_sum = fmaf(a, a, square_sum);
            square_sum = fmaf(b, b, square_sum);
            square_sum = fmaf(c, c, square_sum);
            square_sum = fmaf(d, d, square_sum);
        }
    }

    float total, total_square;
    block_sum_pair(sum, square_sum, total, total_square);
    if (threadIdx.x == 0)
        partial[int64_t(group) * chunks + chunk] =
            make_float2(total, total_square);
}

__global__ void finish_statistics(
    const float2* __restrict__ partial,
    float2* __restrict__ statistics,
    int chunks,
    int64_t elements_per_group,
    float eps) {
    const int group = blockIdx.x;
    float sum = 0.0f;
    float square_sum = 0.0f;

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

    float total, total_square;
    block_sum_pair(sum, square_sum, total, total_square);
    if (threadIdx.x == 0) {
        const float inverse_count = 1.0f / float(elements_per_group);
        const float mean = total * inverse_count;
        float variance = total_square * inverse_count - mean * mean;
        variance = variance < 0.0f ? 0.0f : variance;
        statistics[group] = make_float2(mean, rsqrtf(variance + eps));
    }
}

__global__ void normalize_relu(
    const float4* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float2* __restrict__ statistics,
    float4* __restrict__ output,
    int64_t vectors_per_group,
    int64_t vectors_per_channel,
    int channels_per_group,
    int groups) {
    const int group = blockIdx.y;
    const int first_channel = (group % groups) * channels_per_group;
    const float2 stats = statistics[group];
    const int64_t start = int64_t(blockIdx.x) * OUTPUT_VECTORS;

    #pragma unroll
    for (int i = 0; i < OUTPUT_VECTORS / THREADS; ++i) {
        const int64_t index = start + threadIdx.x + i * THREADS;
        if (index < vectors_per_group) {
            const int channel = first_channel + index / vectors_per_channel;
            const float scale = stats.y * weight[channel];
            const float shift = bias[channel] - stats.x * scale;
            const int64_t offset = int64_t(group) * vectors_per_group + index;
            const float4 v = x[offset];

            output[offset] = make_float4(
                fmaf(relu_value(v.x), scale, shift),
                fmaf(relu_value(v.y), scale, shift),
                fmaf(relu_value(v.z), scale, shift),
                fmaf(relu_value(v.w), 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.scalar_type() == at::kFloat,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");
    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(),
                "Expected contiguous float32 parameters");

    const int64_t batch = x.size(0);
    const int64_t channels = x.size(1);
    const int64_t spatial = x.size(2) * x.size(3) * x.size(4);
    TORCH_CHECK(groups > 0 && channels % groups == 0,
                "Invalid group count");
    TORCH_CHECK(spatial > 0 && spatial % 4 == 0,
                "Spatial size must be a positive multiple of four");
    TORCH_CHECK(weight.numel() == channels && bias.numel() == channels,
                "Invalid parameter sizes");

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

    const int64_t group_count = batch * groups;
    TORCH_CHECK(group_count <= 65535, "Too many groups");
    const int channels_per_group = channels / groups;
    const int64_t elements_per_group = spatial * channels_per_group;
    const int64_t vectors_per_group = elements_per_group / 4;
    const int chunks =
        (vectors_per_group + STATS_VECTORS - 1) / STATS_VECTORS;
    const int output_chunks =
        (vectors_per_group + OUTPUT_VECTORS - 1) / OUTPUT_VECTORS;

    auto partial = torch::empty({group_count, chunks, 2}, x.options());
    auto statistics = torch::empty({group_count, 2}, x.options());
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    partial_statistics<<<dim3(chunks, group_count), THREADS, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        vectors_per_group, chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_statistics<<<group_count, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        reinterpret_cast<float2*>(statistics.data_ptr<float>()),
        chunks, elements_per_group, static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    normalize_relu<<<dim3(output_chunks, group_count), THREADS, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        reinterpret_cast<const float2*>(statistics.data_ptr<float>()),
        reinterpret_cast<float4*>(output.data_ptr<float>()),
        vectors_per_group, spatial / 4, 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)
        spatial = x.shape[2] * x.shape[3] * x.shape[4]
        use_custom = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and spatial > 0
            and spatial % 4 == 0
            and x.shape[0] * self.group_norm.num_groups <= 65535
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
            and not (
                torch.is_grad_enabled()
                and (
                    x.requires_grad
                    or self.group_norm.weight.requires_grad
                    or self.group_norm.bias.requires_grad
                )
            )
        )
        if use_custom:
            return _fused_ops.relu_group_norm_cuda(
                x,
                self.group_norm.weight,
                self.group_norm.bias,
                self.group_norm.num_groups,
                self.group_norm.eps,
            )
        return F.group_norm(
            F.relu(x),
            self.group_norm.num_groups,
            self.group_norm.weight,
            self.group_norm.bias,
            self.group_norm.eps,
        )
```