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

_cpp_source = r"""
#include <torch/extension.h>

torch::Tensor fused_post_conv(
    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 <cmath>

struct Welford {
    float mean;
    float m2;
    int count;
};

__device__ __forceinline__ Welford combine(Welford a, Welford b) {
    if (a.count == 0) return b;
    if (b.count == 0) return a;

    const int count = a.count + b.count;
    const float delta = b.mean - a.mean;
    const float ratio = float(b.count) / float(count);

    Welford result;
    result.mean = a.mean + delta * ratio;
    result.m2 = a.m2 + b.m2
              + delta * delta * (float(a.count) * ratio);
    result.count = count;
    return result;
}

__device__ __forceinline__ Welford warp_reduce(Welford v) {
    for (int offset = 16; offset > 0; offset >>= 1) {
        Welford other;
        other.mean = __shfl_down_sync(0xffffffff, v.mean, offset);
        other.m2 = __shfl_down_sync(0xffffffff, v.m2, offset);
        other.count = __shfl_down_sync(0xffffffff, v.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            v = combine(v, other);
    }
    return v;
}

__global__ void group_statistics(
    const float* __restrict__ x,
    float2* __restrict__ stats,
    int64_t group_elements,
    float eps)
{
    const int group = blockIdx.x;
    const int tid = threadIdx.x;
    const int64_t base = int64_t(group) * group_elements;

    Welford v = {0.0f, 0.0f, 0};
    for (int64_t i = tid; i < group_elements; i += blockDim.x) {
        const float value = x[base + i];
        ++v.count;
        const float delta = value - v.mean;
        v.mean += delta / float(v.count);
        v.m2 += delta * (value - v.mean);
    }

    v = warp_reduce(v);

    __shared__ float warp_mean[8];
    __shared__ float warp_m2[8];
    __shared__ int warp_count[8];

    const int lane = tid & 31;
    const int warp = tid >> 5;
    if (lane == 0) {
        warp_mean[warp] = v.mean;
        warp_m2[warp] = v.m2;
        warp_count[warp] = v.count;
    }
    __syncthreads();

    if (warp == 0) {
        Welford total = {0.0f, 0.0f, 0};
        if (lane < 8) {
            total.mean = warp_mean[lane];
            total.m2 = warp_m2[lane];
            total.count = warp_count[lane];
        }
        total = warp_reduce(total);
        if (lane == 0) {
            const float variance =
                fmaxf(total.m2 / float(total.count), 0.0f);
            stats[group] = make_float2(
                total.mean, rsqrtf(variance + eps));
        }
    }
}

__device__ __forceinline__ float transform(
    float value, float2 stat, float gamma, float beta)
{
    const float normalized = (value - stat.x) * stat.y;
    const float t = tanhf(fmaf(normalized, gamma, beta));

    // tanh lies in [-1, 1], so Hardswish needs no clamping.
    return value + t * (t + 3.0f) * (1.0f / 6.0f);
}

template<int CHANNELS>
__global__ void fused_output(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float2* __restrict__ stats,
    float* __restrict__ output,
    int spatial,
    int groups,
    int channels,
    int64_t pixels)
{
    const int64_t index =
        int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    if (index >= pixels) return;

    const int batch = index / spatial;
    const int position = index % spatial;
    const int c_count = CHANNELS > 0 ? CHANNELS : channels;
    const int channels_per_group = c_count / groups;
    const int64_t base = int64_t(batch) * c_count * spatial + position;
    const float2* batch_stats = stats + int64_t(batch) * groups;

    if constexpr (CHANNELS > 0) {
        float values[CHANNELS];
        float maximum = -CUDART_INF_F;

        #pragma unroll
        for (int c = 0; c < CHANNELS; ++c) {
            const float value = x[base + int64_t(c) * spatial];
            const float2 stat = batch_stats[c / channels_per_group];
            values[c] = transform(value, stat, weight[c], bias[c]);
            maximum = fmaxf(maximum, values[c]);
        }

        float sum = 0.0f;
        #pragma unroll
        for (int c = 0; c < CHANNELS; ++c)
            sum += __expf(values[c] - maximum);

        output[index] = maximum + logf(sum);
    } else {
        // Stable online LogSumExp for other channel counts.
        float maximum = -CUDART_INF_F;
        float sum = 0.0f;
        for (int c = 0; c < channels; ++c) {
            const float value = x[base + int64_t(c) * spatial];
            const float2 stat = batch_stats[c / channels_per_group];
            const float result =
                transform(value, stat, weight[c], bias[c]);
            const float next_maximum = fmaxf(maximum, result);
            sum = sum * __expf(maximum - next_maximum)
                + __expf(result - next_maximum);
            maximum = next_maximum;
        }
        output[index] = maximum + logf(sum);
    }
}

torch::Tensor fused_post_conv(
    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 input");
    TORCH_CHECK(x.dim() == 4 && x.is_contiguous(),
                "Expected contiguous NCHW input");
    TORCH_CHECK(weight.device() == x.device() &&
                bias.device() == x.device() &&
                weight.scalar_type() == at::kFloat &&
                bias.scalar_type() == at::kFloat &&
                weight.is_contiguous() && bias.is_contiguous(),
                "Invalid affine parameters");

    const c10::cuda::CUDAGuard guard(x.device());
    const int64_t batch = x.size(0);
    const int channels = x.size(1);
    const int height = x.size(2);
    const int width = x.size(3);
    const int spatial = height * width;

    TORCH_CHECK(groups > 0 && channels % groups == 0,
                "Invalid group count");
    TORCH_CHECK(weight.numel() == channels && bias.numel() == channels,
                "Incorrect affine parameter size");

    auto output = torch::empty({batch, 1, height, width}, x.options());
    if (batch == 0) return output;

    auto stats = torch::empty({batch * groups, 2}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();
    const int64_t group_elements =
        int64_t(channels / groups) * spatial;

    group_statistics<<<batch * groups, 256, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float2*>(stats.data_ptr<float>()),
        group_elements,
        static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const int64_t pixels = batch * spatial;
    const int blocks = static_cast<int>((pixels + 127) / 128);

    if (channels == 64) {
        fused_output<64><<<blocks, 128, 0, stream>>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            reinterpret_cast<const float2*>(stats.data_ptr<float>()),
            output.data_ptr<float>(),
            spatial, static_cast<int>(groups), channels, pixels);
    } else {
        fused_output<0><<<blocks, 128, 0, stream>>>(
            x.data_ptr<float>(),
            weight.data_ptr<float>(),
            bias.data_ptr<float>(),
            reinterpret_cast<const float2*>(stats.data_ptr<float>()),
            output.data_ptr<float>(),
            spatial, static_cast<int>(groups), channels, pixels);
    }
    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_fused = load_inline(
    name="groupnorm_tanh_hardswish_residual_lse_cuda",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_post_conv"],
    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, eps=1e-5
    ):
        super().__init__()
        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size)
        self.group_norm = nn.GroupNorm(groups, out_channels, eps=eps)
        self.tanh = nn.Tanh()
        self.hard_swish = nn.Hardswish()

    def forward(self, x):
        x_conv = self.conv(x)

        if (
            x_conv.is_cuda
            and x_conv.dtype == torch.float32
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
        ):
            return _fused.fused_post_conv(
                x_conv.contiguous(),
                self.group_norm.weight.contiguous(),
                self.group_norm.bias.contiguous(),
                self.group_norm.num_groups,
                self.group_norm.eps,
            )

        x_norm = self.group_norm(x_conv)
        x_res = x_conv + self.hard_swish(self.tanh(x_norm))
        return torch.logsumexp(x_res, dim=1, keepdim=True)
```