```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_groupnorm_lse_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.h>
#include <cuda_runtime.h>
#include <cmath>

__device__ __forceinline__ void merge_welford(
    float& mean, float& m2, int& count,
    float other_mean, float other_m2, int other_count)
{
    const int total = count + other_count;
    if (other_count && total) {
        const float delta = other_mean - mean;
        const float ratio = float(other_count) / float(total);
        m2 += other_m2 + delta * delta * float(count) * ratio;
        mean += delta * ratio;
    }
    count = total;
}

__device__ __forceinline__ void warp_welford(
    float& mean, float& m2, int& count)
{
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        const float other_mean =
            __shfl_down_sync(0xffffffff, mean, offset);
        const float other_m2 =
            __shfl_down_sync(0xffffffff, m2, offset);
        const int other_count =
            __shfl_down_sync(0xffffffff, count, offset);

        if ((threadIdx.x & 31) + offset < 32) {
            merge_welford(
                mean, m2, count, other_mean, other_m2, other_count);
        }
    }
}

__global__ void group_statistics(
    const float* __restrict__ x,
    float2* __restrict__ statistics,
    int64_t group_elements,
    float eps)
{
    __shared__ float warp_means[8];
    __shared__ float warp_m2[8];
    __shared__ int warp_counts[8];

    const int tid = threadIdx.x;
    const int lane = tid & 31;
    const int warp = tid >> 5;
    const int64_t base = int64_t(blockIdx.x) * group_elements;

    float mean = 0.0f;
    float m2 = 0.0f;
    int count = 0;

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

    warp_welford(mean, m2, count);

    if (lane == 0) {
        warp_means[warp] = mean;
        warp_m2[warp] = m2;
        warp_counts[warp] = count;
    }
    __syncthreads();

    if (warp == 0) {
        mean = lane < 8 ? warp_means[lane] : 0.0f;
        m2 = lane < 8 ? warp_m2[lane] : 0.0f;
        count = lane < 8 ? warp_counts[lane] : 0;

        warp_welford(mean, m2, count);

        if (lane == 0) {
            const float variance =
                fmaxf(m2 / float(group_elements), 0.0f);
            statistics[blockIdx.x] =
                make_float2(mean, rsqrtf(variance + eps));
        }
    }
}

__global__ void fused_activation_residual_lse64(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float2* __restrict__ statistics,
    float* __restrict__ output,
    int64_t spatial_size,
    int groups)
{
    const int64_t pixel =
        int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
    if (pixel >= spatial_size)
        return;

    const int batch = blockIdx.y;
    const int channels_per_group = 64 / groups;
    const int64_t base = int64_t(batch) * 64 * spatial_size + pixel;

    float values[64];
    float maximum = -CUDART_INF_F;

    #pragma unroll
    for (int channel = 0; channel < 64; ++channel) {
        const float residual = x[base + int64_t(channel) * spatial_size];
        const int group = channel / channels_per_group;
        const float2 stats = statistics[batch * groups + group];

        const float normalized =
            ((residual - stats.x) * stats.y) * weight[channel]
            + bias[channel];
        const float activated = tanhf(normalized);

        // Tanh lies in [-1, 1], so the Hardswish clamp is unnecessary.
        const float value =
            residual + activated * (activated + 3.0f) / 6.0f;

        values[channel] = value;
        maximum = fmaxf(maximum, value);
    }

    float sum = 0.0f;
    #pragma unroll
    for (int channel = 0; channel < 64; ++channel) {
        sum += __expf(values[channel] - maximum);
    }

    output[int64_t(batch) * spatial_size + pixel] =
        maximum + logf(sum);
}

torch::Tensor fused_groupnorm_lse_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() == 4 && x.size(1) == 64 && x.is_contiguous(),
                "Expected contiguous NCHW input with 64 channels");
    TORCH_CHECK(groups > 0 && 64 % groups == 0,
                "Invalid group count");
    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()
                && weight.numel() == 64
                && bias.numel() == 64,
                "Invalid affine parameters");

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

    const int64_t batch = x.size(0);
    const int64_t height = x.size(2);
    const int64_t width = x.size(3);
    const int64_t spatial_size = height * width;

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

    auto statistics = torch::empty({batch, groups, 2}, x.options());
    const auto stream = at::cuda::getCurrentCUDAStream();

    group_statistics<<<batch * groups, 256, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float2*>(statistics.data_ptr<float>()),
        (64 / groups) * spatial_size,
        static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    constexpr int threads = 128;
    const dim3 grid(
        static_cast<unsigned int>((spatial_size + threads - 1) / threads),
        static_cast<unsigned int>(batch));

    fused_activation_residual_lse64<<<grid, threads, 0, stream>>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        reinterpret_cast<const float2*>(statistics.data_ptr<float>()),
        output.data_ptr<float>(),
        spatial_size,
        static_cast<int>(groups));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_fused_ops = load_inline(
    name="fused_groupnorm_tanh_hardswish_residual_lse64_v1",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["fused_groupnorm_lse_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    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 x_conv.shape[1] == 64
            and x_conv.is_contiguous()
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
        ):
            return _fused_ops.fused_groupnorm_lse_cuda(
                x_conv,
                self.group_norm.weight,
                self.group_norm.bias,
                self.group_norm.num_groups,
                self.group_norm.eps,
            )

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