```python
import torch
import torch.nn as nn
import torch.nn.functional as F
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 <stdint.h>

constexpr int THREADS = 256;
constexpr int TILE = 4096;

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

__device__ __forceinline__ void reduce_pair(
    float& sum, float& square_sum
) {
    __shared__ float warp_sums[8];
    __shared__ float warp_squares[8];

    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) {
        warp_sums[warp] = sum;
        warp_squares[warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = lane < 8 ? warp_sums[lane] : 0.0f;
        square_sum = lane < 8 ? warp_squares[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
            );
        }
    }
}

__global__ void relu_partial_stats(
    const float* __restrict__ x,
    float* __restrict__ partial,
    int64_t group_size,
    int tiles
) {
    const int group = blockIdx.x / tiles;
    const int tile = blockIdx.x % tiles;
    const int64_t start = int64_t(tile) * TILE;
    const int64_t base = int64_t(group) * group_size;

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

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

    reduce_pair(sum, square_sum);

    if (threadIdx.x == 0) {
        const int64_t index = int64_t(blockIdx.x) * 2;
        partial[index] = sum;
        partial[index + 1] = square_sum;
    }
}

__global__ void finish_stats(
    const float* __restrict__ partial,
    float* __restrict__ stats,
    int64_t group_size,
    int tiles,
    float eps
) {
    const int group = blockIdx.x;
    float sum = 0.0f;
    float square_sum = 0.0f;

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

    reduce_pair(sum, square_sum);

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

__global__ void relu_normalize_affine(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float* __restrict__ stats,
    float* __restrict__ output,
    int64_t group_size,
    int64_t spatial_size,
    int channels_per_group,
    int groups,
    int tiles
) {
    const int group = blockIdx.x / tiles;
    const int tile = blockIdx.x % tiles;
    const int64_t start = int64_t(tile) * TILE;
    const int64_t base = int64_t(group) * group_size;
    const int first_channel = (group % groups) * channels_per_group;

    const float mean = stats[group * 2];
    const float inverse_std = stats[group * 2 + 1];

    #pragma unroll
    for (int k = 0; k < TILE / THREADS; ++k) {
        const int64_t offset = start + threadIdx.x + k * THREADS;
        if (offset < group_size) {
            const int channel = first_channel + offset / spatial_size;
            const float value = relu_value(x[base + offset]);
            const float normalized = (value - mean) * inverse_std;
            output[base + offset] =
                normalized * weight[channel] + bias[channel];
        }
    }
}

torch::Tensor fused_relu_group_norm(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps
) {
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat);
    TORCH_CHECK(x.is_contiguous() && x.dim() == 5);
    TORCH_CHECK(weight.device() == x.device());
    TORCH_CHECK(bias.device() == x.device());
    TORCH_CHECK(weight.scalar_type() == at::kFloat);
    TORCH_CHECK(bias.scalar_type() == at::kFloat);
    TORCH_CHECK(weight.is_contiguous() && bias.is_contiguous());

    const c10::cuda::CUDAGuard device_guard(x.device());

    const int64_t channels = x.size(1);
    TORCH_CHECK(groups > 0 && channels % groups == 0);
    TORCH_CHECK(weight.numel() == channels && bias.numel() == channels);

    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 int64_t group_size = (channels / groups) * spatial_size;
    const int64_t total_groups = x.size(0) * groups;
    const int tiles = (group_size + TILE - 1) / TILE;

    auto partial = torch::empty({total_groups, tiles, 2}, x.options());
    auto stats = torch::empty({total_groups, 2}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    relu_partial_stats<<<total_groups * tiles, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        partial.data_ptr<float>(),
        group_size,
        tiles
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_stats<<<total_groups, THREADS, 0, stream>>>(
        partial.data_ptr<float>(),
        stats.data_ptr<float>(),
        group_size,
        tiles,
        static_cast<float>(eps)
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    relu_normalize_affine<<<total_groups * tiles, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        stats.data_ptr<float>(),
        output.data_ptr<float>(),
        group_size,
        spatial_size,
        channels / groups,
        groups,
        tiles
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_extension = load_inline(
    name="fused_relu_group_norm_3d_v1",
    cpp_sources="""
        torch::Tensor fused_relu_group_norm(
            torch::Tensor x,
            torch::Tensor weight,
            torch::Tensor bias,
            int64_t groups,
            double eps
        );
    """,
    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
        )
        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)
        norm = self.group_norm

        needs_grad = torch.is_grad_enabled() and (
            x.requires_grad
            or norm.weight.requires_grad
            or norm.bias.requires_grad
        )

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and norm.weight.dtype == torch.float32
            and norm.bias.dtype == torch.float32
            and not needs_grad
        ):
            return _extension.fused_relu_group_norm(
                x,
                norm.weight,
                norm.bias,
                norm.num_groups,
                norm.eps,
            )

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