```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 pack_input_cuda(torch::Tensor input);

torch::Tensor fused_post_conv_cuda(
    torch::Tensor input,
    torch::Tensor conv_bias,
    torch::Tensor weight,
    torch::Tensor bias,
    double eps);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAStream.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <math_constants.h>

constexpr int CHANNELS = 64;
constexpr int GROUPS = 16;
constexpr int CHUNK = 512;

__device__ __forceinline__ float warp_sum(float x) {
    #pragma unroll
    for (int d = 16; d; d >>= 1)
        x += __shfl_down_sync(0xffffffffu, x, d);
    return x;
}

__device__ __forceinline__ float fast_tanh(float x) {
#if __CUDA_ARCH__ >= 750
    float y;
    asm("tanh.approx.f32 %0, %1;" : "=f"(y) : "f"(x));
    return y;
#else
    return tanhf(x);
#endif
}

__global__ void pack_eight_channels(
    const float* __restrict__ input,
    float* __restrict__ output,
    int spatial
) {
    __shared__ float tile[8][257];

    const int tid = threadIdx.x;
    const int start = blockIdx.x * 256;
    const int pixel = start + tid;
    const int batch = blockIdx.y;
    const int64_t input_base = static_cast<int64_t>(batch) * 8 * spatial;

    #pragma unroll
    for (int c = 0; c < 8; ++c) {
        tile[c][tid] = pixel < spatial
            ? input[input_base + static_cast<int64_t>(c) * spatial + pixel]
            : 0.0f;
    }
    __syncthreads();

    #pragma unroll
    for (int part = 0; part < 2; ++part) {
        const int vector = tid + part * 256;
        const int p = vector >> 1;
        const int c = (vector & 1) * 4;
        if (start + p < spatial) {
            const float4 value = make_float4(
                tile[c][p], tile[c + 1][p],
                tile[c + 2][p], tile[c + 3][p]
            );
            reinterpret_cast<float4*>(output)[
                (static_cast<int64_t>(batch) * spatial + start + p) * 2
                + (vector & 1)
            ] = value;
        }
    }
}

torch::Tensor pack_input_cuda(torch::Tensor input) {
    TORCH_CHECK(
        input.is_cuda() && input.scalar_type() == torch::kFloat32 &&
        input.dim() == 4 && input.size(1) == 8 && input.is_contiguous(),
        "Expected contiguous NCHW CUDA float32 input with eight channels"
    );

    const c10::cuda::CUDAGuard guard(input.device());
    auto output = torch::empty(
        input.sizes(),
        input.options().memory_format(at::MemoryFormat::ChannelsLast)
    );

    const int spatial = static_cast<int>(input.size(2) * input.size(3));
    const int batch = static_cast<int>(input.size(0));

    if (batch && spatial) {
        const auto stream =
            c10::cuda::getCurrentCUDAStream(input.get_device()).stream();

        pack_eight_channels<<<
            dim3((spatial + 255) / 256, batch), 256, 0, stream
        >>>(
            input.data_ptr<float>(), output.data_ptr<float>(), spatial
        );
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }
    return output;
}

__global__ void partial_moments(
    const float* __restrict__ input,
    const float* __restrict__ conv_bias,
    float2* __restrict__ partial,
    int spatial,
    int chunks
) {
    const int tid = threadIdx.x;
    const int group = tid & 15;
    const int spatial_lane = tid >> 4;
    const int batch = blockIdx.y;
    const int start = blockIdx.x * CHUNK;

    const float4 b = reinterpret_cast<const float4*>(conv_bias)[group];

    float sx = 0.0f, sy = 0.0f, sz = 0.0f, sw = 0.0f;
    float qx = 0.0f, qy = 0.0f, qz = 0.0f, qw = 0.0f;

    #pragma unroll 4
    for (int offset = 0; offset < CHUNK; offset += 16) {
        const int position = start + spatial_lane + offset;
        if (position < spatial) {
            float4 v = reinterpret_cast<const float4*>(input)[
                (static_cast<int64_t>(batch) * spatial + position) * GROUPS
                + group
            ];
            v.x += b.x;
            v.y += b.y;
            v.z += b.z;
            v.w += b.w;

            sx += v.x;
            sy += v.y;
            sz += v.z;
            sw += v.w;

            qx = fmaf(v.x, v.x, qx);
            qy = fmaf(v.y, v.y, qy);
            qz = fmaf(v.z, v.z, qz);
            qw = fmaf(v.w, v.w, qw);
        }
    }

    float sum = (sx + sy) + (sz + sw);
    float square_sum = (qx + qy) + (qz + qw);

    sum += __shfl_xor_sync(0xffffffffu, sum, 16);
    square_sum += __shfl_xor_sync(0xffffffffu, square_sum, 16);

    __shared__ float shared_sum[8][16];
    __shared__ float shared_square[8][16];

    if ((tid & 31) < 16) {
        shared_sum[tid >> 5][group] = sum;
        shared_square[tid >> 5][group] = square_sum;
    }
    __syncthreads();

    if (tid < GROUPS) {
        float total = 0.0f;
        float total_square = 0.0f;

        #pragma unroll
        for (int w = 0; w < 8; ++w) {
            total += shared_sum[w][tid];
            total_square += shared_square[w][tid];
        }

        partial[
            (static_cast<int64_t>(batch) * chunks + blockIdx.x) * GROUPS
            + tid
        ] = make_float2(total, total_square);
    }
}

__global__ void finish_moments(
    const float2* __restrict__ partial,
    float2* __restrict__ moments,
    int spatial,
    int chunks,
    float eps
) {
    const int row = blockIdx.x;
    const int batch = row / GROUPS;
    const int group = row % GROUPS;

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

    for (int chunk = threadIdx.x; chunk < chunks; chunk += 32) {
        const float2 v = partial[
            (static_cast<int64_t>(batch) * chunks + chunk) * GROUPS + group
        ];
        sum += v.x;
        square_sum += v.y;
    }

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

    if (threadIdx.x == 0) {
        const float inverse_count = 1.0f / (4.0f * spatial);
        const float mean = sum * inverse_count;
        const float variance =
            fmaxf(fmaf(-mean, mean, square_sum * inverse_count), 0.0f);
        moments[row] = make_float2(mean, rsqrtf(variance + eps));
    }
}

__device__ __forceinline__ float activate(
    float x, float cb, float weight, float bias, float2 stats
) {
    const float original = x + cb;
    const float normalized =
        fmaf((original - stats.x) * stats.y, weight, bias);
    const float t = fast_tanh(normalized);
    const float h = (t * (t + 3.0f)) * (1.0f / 6.0f);
    return original + h;
}

__global__ void fused_output(
    const float* __restrict__ input,
    const float* __restrict__ conv_bias,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    const float2* __restrict__ moments,
    float* __restrict__ output,
    int spatial
) {
    const int pixel = blockIdx.x * 32 + (threadIdx.x >> 3);
    const int lane = threadIdx.x & 7;
    const int batch = blockIdx.y;

    if (pixel >= spatial)
        return;

    const unsigned mask = __activemask();
    const int64_t base =
        (static_cast<int64_t>(batch) * spatial + pixel) * GROUPS;

    float values[8];

    #pragma unroll
    for (int half = 0; half < 2; ++half) {
        // Each load covers eight consecutive float4 vectors.
        const int group = lane + half * 8;
        const float4 x =
            reinterpret_cast<const float4*>(input)[base + group];
        const float4 cb =
            reinterpret_cast<const float4*>(conv_bias)[group];
        const float4 w =
            reinterpret_cast<const float4*>(weight)[group];
        const float4 b =
            reinterpret_cast<const float4*>(bias)[group];
        const float2 stats = moments[batch * GROUPS + group];

        values[half * 4] =
            activate(x.x, cb.x, w.x, b.x, stats);
        values[half * 4 + 1] =
            activate(x.y, cb.y, w.y, b.y, stats);
        values[half * 4 + 2] =
            activate(x.z, cb.z, w.z, b.z, stats);
        values[half * 4 + 3] =
            activate(x.w, cb.w, w.w, b.w, stats);
    }

    float maximum = -CUDART_INF_F;
    #pragma unroll
    for (int i = 0; i < 8; ++i)
        maximum = fmaxf(maximum, values[i]);

    #pragma unroll
    for (int d = 4; d; d >>= 1)
        maximum = fmaxf(
            maximum, __shfl_xor_sync(mask, maximum, d, 8)
        );

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

    #pragma unroll
    for (int d = 4; d; d >>= 1)
        sum += __shfl_xor_sync(mask, sum, d, 8);

    if (lane == 0)
        output[static_cast<int64_t>(batch) * spatial + pixel] =
            maximum + __logf(sum);
}

torch::Tensor fused_post_conv_cuda(
    torch::Tensor input,
    torch::Tensor conv_bias,
    torch::Tensor weight,
    torch::Tensor bias,
    double eps
) {
    TORCH_CHECK(
        input.is_cuda() && input.scalar_type() == torch::kFloat32 &&
        input.dim() == 4 && input.size(1) == CHANNELS &&
        input.is_contiguous(at::MemoryFormat::ChannelsLast),
        "Expected channels-last CUDA float32 input with 64 channels"
    );

    for (const auto& p : {conv_bias, weight, bias}) {
        TORCH_CHECK(
            p.device() == input.device() &&
            p.scalar_type() == torch::kFloat32 &&
            p.is_contiguous() && p.numel() == CHANNELS,
            "Expected 64 contiguous float32 parameters on the input device"
        );
    }

    const c10::cuda::CUDAGuard guard(input.device());
    const auto stream =
        c10::cuda::getCurrentCUDAStream(input.get_device()).stream();

    const int batch = static_cast<int>(input.size(0));
    const int spatial = static_cast<int>(input.size(2) * input.size(3));

    auto output = torch::empty(
        {input.size(0), 1, input.size(2), input.size(3)}, input.options()
    );
    if (!batch || !spatial)
        return output;

    const int chunks = (spatial + CHUNK - 1) / CHUNK;
    auto partial = torch::empty({batch, chunks, GROUPS, 2}, input.options());
    auto moments = torch::empty({batch, GROUPS, 2}, input.options());

    partial_moments<<<dim3(chunks, batch), 256, 0, stream>>>(
        input.data_ptr<float>(),
        conv_bias.data_ptr<float>(),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        spatial, chunks
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_moments<<<batch * GROUPS, 32, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        reinterpret_cast<float2*>(moments.data_ptr<float>()),
        spatial, chunks, static_cast<float>(eps)
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    fused_output<<<dim3((spatial + 31) / 32, batch), 256, 0, stream>>>(
        input.data_ptr<float>(),
        conv_bias.data_ptr<float>(),
        weight.data_ptr<float>(),
        bias.data_ptr<float>(),
        reinterpret_cast<const float2*>(moments.data_ptr<float>()),
        output.data_ptr<float>(),
        spatial
    );
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

_ops = load_inline(
    name="fused_nhwc_bias_gn_tanh_hswish_residual_lse_v5",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["pack_input_cuda", "fused_post_conv_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.conv.to(memory_format=torch.channels_last)
        self.group_norm = nn.GroupNorm(groups, out_channels, eps=eps)
        self.tanh = nn.Tanh()
        self.hard_swish = nn.Hardswish()

    def forward(self, x):
        if (
            x.is_cuda
            and x.dtype == torch.float32
            and self.conv.out_channels == 64
            and self.group_norm.num_groups == 16
            and self.conv.weight.dtype == torch.float32
            and self.conv.bias is not None
            and self.conv.bias.dtype == torch.float32
            and self.group_norm.weight.dtype == torch.float32
            and self.group_norm.bias.dtype == torch.float32
            and not torch.is_autocast_enabled()
        ):
            if not x.is_contiguous(memory_format=torch.channels_last):
                if x.shape[1] == 8 and x.is_contiguous():
                    x = _ops.pack_input_cuda(x)
                else:
                    x = x.contiguous(memory_format=torch.channels_last)

            weight = self.conv.weight.contiguous(
                memory_format=torch.channels_last
            )
            x_conv = F.conv2d(
                x,
                weight,
                None,
                self.conv.stride,
                self.conv.padding,
                self.conv.dilation,
                self.conv.groups,
            )
            return _ops.fused_post_conv_cuda(
                x_conv,
                self.conv.bias,
                self.group_norm.weight,
                self.group_norm.bias,
                self.group_norm.eps,
            )

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