```python
import torch
import torch.nn as nn
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>

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

__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;

    #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<bool VECTORIZE>
__global__ void relu_partial_stats(
    const float *__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 start = int64_t(tile) * TILE;
    const float *input = x + int64_t(group) * group_size;
    float sum = 0.0f;
    float square_sum = 0.0f;

    if constexpr (VECTORIZE) {
        const float4 *input4 = reinterpret_cast<const float4 *>(input);
        #pragma unroll
        for (int j = 0; j < TILE / (4 * THREADS); ++j) {
            const int64_t i = start / 4 + threadIdx.x + j * THREADS;
            if (i < group_size / 4) {
                float4 v = input4[i];
                v.x = fmaxf(v.x, 0.0f);
                v.y = fmaxf(v.y, 0.0f);
                v.z = fmaxf(v.z, 0.0f);
                v.w = fmaxf(v.w, 0.0f);
                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;
            }
        }
    } else {
        #pragma unroll
        for (int j = 0; j < TILE / THREADS; ++j) {
            const int64_t i = start + threadIdx.x + j * THREADS;
            if (i < group_size) {
                const float v = fmaxf(input[i], 0.0f);
                sum += v;
                square_sum += v * v;
            }
        }
    }

    block_sum(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,
    int tiles,
    int64_t group_size,
    float eps) {
    __shared__ float shared_sum[8];
    __shared__ float shared_square_sum[8];

    float sum = 0.0f;
    float square_sum = 0.0f;
    const int64_t base = int64_t(blockIdx.x) * tiles;

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

    block_sum(sum, square_sum, shared_sum, shared_square_sum);
    if (threadIdx.x == 0) {
        const float mean = sum / float(group_size);
        const float variance =
            fmaxf(square_sum / float(group_size) - mean * mean, 0.0f);
        stats[blockIdx.x] = make_float2(mean, rsqrtf(variance + eps));
    }
}

template<bool VECTORIZE>
__global__ void relu_normalize(
    const float *__restrict__ x,
    const float *__restrict__ weight,
    const float *__restrict__ bias,
    const float2 *__restrict__ stats,
    float *__restrict__ output,
    int64_t group_size,
    int64_t spatial,
    int groups,
    int channels_per_group) {
    const int group = blockIdx.y;
    const int64_t start = int64_t(blockIdx.x) * TILE;
    const int64_t base = int64_t(group) * group_size;
    const int first_channel = (group % groups) * channels_per_group;
    const float2 stat = stats[group];

    if constexpr (VECTORIZE) {
        const float4 *input4 =
            reinterpret_cast<const float4 *>(x + base);
        float4 *output4 = reinterpret_cast<float4 *>(output + base);

        #pragma unroll
        for (int j = 0; j < TILE / (4 * THREADS); ++j) {
            const int64_t i = start / 4 + threadIdx.x + j * THREADS;
            if (i < group_size / 4) {
                const int channel = first_channel + (4 * i) / spatial;
                const float scale = stat.y * weight[channel];
                const float shift = bias[channel] - stat.x * scale;
                float4 v = input4[i];
                v.x = fmaf(fmaxf(v.x, 0.0f), scale, shift);
                v.y = fmaf(fmaxf(v.y, 0.0f), scale, shift);
                v.z = fmaf(fmaxf(v.z, 0.0f), scale, shift);
                v.w = fmaf(fmaxf(v.w, 0.0f), scale, shift);
                output4[i] = v;
            }
        }
    } else {
        #pragma unroll
        for (int j = 0; j < TILE / THREADS; ++j) {
            const int64_t i = start + threadIdx.x + j * THREADS;
            if (i < group_size) {
                const int channel = first_channel + i / spatial;
                const float scale = stat.y * weight[channel];
                const float shift = bias[channel] - stat.x * scale;
                output[base + i] =
                    fmaf(fmaxf(x[base + i], 0.0f), scale, shift);
            }
        }
    }
}

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.dim() == 5 && x.is_contiguous());
    TORCH_CHECK(weight.is_cuda() && bias.is_cuda());
    TORCH_CHECK(weight.device() == x.device() && bias.device() == x.device());
    TORCH_CHECK(weight.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat);
    TORCH_CHECK(weight.is_contiguous() && bias.is_contiguous());
    TORCH_CHECK(groups > 0 && x.size(1) % groups == 0);
    TORCH_CHECK(weight.numel() == x.size(1) && bias.numel() == x.size(1));

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

    const int64_t spatial = 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;
    const int64_t total_groups = x.size(0) * groups;
    TORCH_CHECK(total_groups <= 65535);

    const int tiles = (group_size + TILE - 1) / TILE;
    auto partial = torch::empty({total_groups, int64_t(tiles), 2}, x.options());
    auto stats = torch::empty({total_groups, 2}, x.options());

    const dim3 grid(tiles, total_groups);
    const auto stream = at::cuda::getCurrentCUDAStream();
    auto partial_ptr = reinterpret_cast<float2 *>(partial.data_ptr<float>());
    auto stats_ptr = reinterpret_cast<float2 *>(stats.data_ptr<float>());

    if (group_size % 4 == 0) {
        relu_partial_stats<true><<<grid, THREADS, 0, stream>>>(
            x.data_ptr<float>(), partial_ptr, group_size, tiles);
    } else {
        relu_partial_stats<false><<<grid, THREADS, 0, stream>>>(
            x.data_ptr<float>(), partial_ptr, group_size, tiles);
    }

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

    if (spatial % 4 == 0) {
        relu_normalize<true><<<grid, THREADS, 0, stream>>>(
            x.data_ptr<float>(), weight.data_ptr<float>(), bias.data_ptr<float>(),
            stats_ptr, output.data_ptr<float>(), group_size, spatial,
            int(groups), channels_per_group);
    } else {
        relu_normalize<false><<<grid, THREADS, 0, stream>>>(
            x.data_ptr<float>(), weight.data_ptr<float>(), bias.data_ptr<float>(),
            stats_ptr, output.data_ptr<float>(), group_size, spatial,
            int(groups), channels_per_group);
    }

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused_ops = load_inline(
    name="fused_relu_group_norm_3d_v1",
    cpp_sources="""
    #include <torch/extension.h>
    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", "-std=c++17"],
    extra_cuda_cflags=["-O3", "-std=c++17"],
    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)
        gn = self.group_norm

        if (
            x.is_cuda
            and x.dtype == torch.float32
            and gn.weight.dtype == torch.float32
            and gn.bias.dtype == torch.float32
            and x.shape[0] * gn.num_groups <= 65535
        ):
            return _fused_ops.fused_relu_group_norm(
                x.contiguous(),
                gn.weight.contiguous(),
                gn.bias.contiguous(),
                gn.num_groups,
                gn.eps,
            )

        return gn(self.relu(x))
```