```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>
#include <algorithm>

constexpr int THREADS = 256;
constexpr int VECTORS_PER_CHUNK = 1024;

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

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

    sum = warp_sum(sum);
    square_sum = warp_sum(square_sum);

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

    if (warp == 0) {
        sum = lane < THREADS / 32 ? shared_sum[lane] : 0.0f;
        square_sum =
            lane < THREADS / 32 ? shared_square_sum[lane] : 0.0f;
        sum = warp_sum(sum);
        square_sum = warp_sum(square_sum);
    }
}

__global__ void relu_partial_stats(
    const float4* __restrict__ x,
    float2* __restrict__ partial,
    int64_t vectors_per_group,
    int chunks) {
    __shared__ float shared_sum[THREADS / 32];
    __shared__ float shared_square_sum[THREADS / 32];

    const int group = blockIdx.x;
    const int chunk = blockIdx.y;
    const int64_t start =
        static_cast<int64_t>(chunk) * VECTORS_PER_CHUNK;
    const int64_t base =
        static_cast<int64_t>(group) * vectors_per_group;

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

    #pragma unroll
    for (int i = 0; i < VECTORS_PER_CHUNK / THREADS; ++i) {
        const int64_t index = start + threadIdx.x + i * THREADS;
        if (index < vectors_per_group) {
            float4 v = x[base + index];
            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;

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

    block_sum(sum, square_sum, shared_sum, shared_square_sum);
    if (threadIdx.x == 0)
        partial[static_cast<int64_t>(group) * chunks + chunk] =
            make_float2(sum, square_sum);
}

__global__ void finish_stats(
    const float2* __restrict__ partial,
    float2* __restrict__ stats,
    int chunks,
    int64_t elements_per_group,
    float eps) {
    __shared__ float shared_sum[THREADS / 32];
    __shared__ float shared_square_sum[THREADS / 32];

    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 p =
            partial[static_cast<int64_t>(group) * chunks + i];
        sum += p.x;
        square_sum += p.y;
    }

    block_sum(sum, square_sum, shared_sum, shared_square_sum);
    if (threadIdx.x == 0) {
        const float reciprocal_count =
            1.0f / static_cast<float>(elements_per_group);
        const float mean = sum * reciprocal_count;
        const float variance =
            fmaxf(square_sum * reciprocal_count - mean * mean, 0.0f);
        stats[group] = make_float2(mean, rsqrtf(variance + eps));
    }
}

__global__ void relu_normalize(
    const float4* __restrict__ x,
    float4* __restrict__ out,
    const float2* __restrict__ stats,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    int64_t vector_count,
    int64_t vectors_per_group,
    int64_t vectors_per_channel,
    int channels) {
    const int64_t stride =
        static_cast<int64_t>(blockDim.x) * gridDim.x;

    for (int64_t index =
             static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
         index < vector_count;
         index += stride) {
        const int64_t group = index / vectors_per_group;
        const int channel = (index / vectors_per_channel) % channels;

        const float2 st = stats[group];
        const float scale = st.y * weight[channel];
        const float shift = bias[channel];

        float4 v = x[index];
        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;

        v.x = (v.x - st.x) * scale + shift;
        v.y = (v.y - st.x) * scale + shift;
        v.z = (v.z - st.x) * scale + shift;
        v.w = (v.w - st.x) * scale + shift;
        out[index] = v;
    }
}

torch::Tensor relu_group_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    int64_t groups,
    double eps) {
    const c10::cuda::CUDAGuard device_guard(x.device());

    const int64_t channels = x.size(1);
    const int64_t spatial = x.size(2) * x.size(3) * x.size(4);
    const int64_t elements_per_group = channels / groups * spatial;
    const int64_t vectors_per_group = elements_per_group / 4;
    const int64_t total_groups = x.size(0) * groups;
    const int chunks =
        (vectors_per_group + VECTORS_PER_CHUNK - 1) / VECTORS_PER_CHUNK;

    auto output = torch::empty(x.sizes(), x.options());
    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();

    relu_partial_stats<<<dim3(total_groups, chunks), 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_stats<<<total_groups, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        reinterpret_cast<float2*>(stats.data_ptr<float>()),
        chunks,
        elements_per_group,
        static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const int64_t vector_count = x.numel() / 4;
    const int blocks = static_cast<int>(
        std::min<int64_t>((vector_count + THREADS - 1) / THREADS, 8192));

    relu_normalize<<<blocks, 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>()),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        vector_count,
        vectors_per_group,
        spatial / 4,
        static_cast<int>(channels));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_fused_ops = load_inline(
    name="fused_relu_groupnorm_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(groups, out_channels)

    def forward(self, x):
        x = self.conv_transpose(x)
        norm = self.group_norm
        spatial = x.shape[2] * x.shape[3] * x.shape[4]
        group_elements = (x.shape[1] // norm.num_groups) * spatial

        use_custom = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.is_contiguous()
            and x.numel() > 0
            and spatial % 4 == 0
            and (group_elements + 4095) // 4096 <= 65535
            and norm.weight.dtype == torch.float32
            and norm.bias.dtype == torch.float32
            and norm.weight.is_contiguous()
            and norm.bias.is_contiguous()
            and not (
                torch.is_grad_enabled()
                and (
                    x.requires_grad
                    or norm.weight.requires_grad
                    or norm.bias.requires_grad
                )
            )
        )

        if use_custom:
            return _fused_ops.relu_group_norm_cuda(
                x, norm.weight, norm.bias, norm.num_groups, norm.eps
            )

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