```python
# Copyright (c) Meta Platforms, Inc. and affiliates.
# All rights reserved.
#
# This source code is licensed under the license found in the
# LICENSE file in the root directory of this source tree.

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

# Custom CUDA kernel that fuses: add_value + Swish + Tanh + GELU + Hardtanh
# into a single element-wise pass over the matmul output.
fused_ops_source = """
#include <torch/extension.h>
#include <cuda_runtime.h>

__device__ __forceinline__ float apply_ops(float v) {
    // Swish: x * sigmoid(x)
    float sig = 1.0f / (1.0f + expf(-v));
    float s = v * sig;
    // Tanh
    float t = tanhf(s);
    // GELU (exact, erf-based, matches PyTorch default)
    float g = 0.5f * t * (1.0f + erff(t * 0.70710678118654752440f));
    // Hardtanh clamp to [-1, 1]
    return fminf(fmaxf(g, -1.0f), 1.0f);
}

__global__ void fused_ops_vec4_kernel(const float4* __restrict__ x,
                                      const float4* __restrict__ add_val,
                                      float4* __restrict__ out,
                                      int cols4, int rows) {
    int col4 = blockIdx.x * blockDim.x + threadIdx.x;
    if (col4 >= cols4) return;
    for (int row = blockIdx.y; row < rows; row += gridDim.y) {
        long idx = (long)row * cols4 + col4;
        float4 xv = x[idx];
        float4 av = add_val[col4];
        float4 res;
        res.x = apply_ops(xv.x + av.x);
        res.y = apply_ops(xv.y + av.y);
        res.z = apply_ops(xv.z + av.z);
        res.w = apply_ops(xv.w + av.w);
        out[idx] = res;
    }
}

__global__ void fused_ops_scalar_kernel(const float* __restrict__ x,
                                        const float* __restrict__ add_val,
                                        float* __restrict__ out,
                                        int cols, int rows) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (col >= cols) return;
    for (int row = blockIdx.y; row < rows; row += gridDim.y) {
        long idx = (long)row * cols + col;
        out[idx] = apply_ops(x[idx] + add_val[col]);
    }
}

torch::Tensor fused_ops_cuda(torch::Tensor x, torch::Tensor add_val) {
    TORCH_CHECK(x.is_cuda() && add_val.is_cuda(), "Inputs must be CUDA tensors");
    TORCH_CHECK(x.scalar_type() == torch::kFloat32, "x must be float32");
    TORCH_CHECK(add_val.scalar_type() == torch::kFloat32, "add_val must be float32");

    auto x_c = x.contiguous();
    auto add_c = add_val.contiguous();
    auto out = torch::empty_like(x_c);

    int cols = x_c.size(-1);
    TORCH_CHECK(add_c.numel() == cols, "add_val must match the last dimension of x");
    long total = x_c.numel();
    int rows = (int)(total / cols);

    const int block_size = 256;
    int grid_y = rows < 65535 ? rows : 65535;

    if (cols % 4 == 0) {
        int cols4 = cols / 4;
        dim3 grid((cols4 + block_size - 1) / block_size, grid_y);
        fused_ops_vec4_kernel<<<grid, block_size>>>(
            reinterpret_cast<const float4*>(x_c.data_ptr<float>()),
            reinterpret_cast<const float4*>(add_c.data_ptr<float>()),
            reinterpret_cast<float4*>(out.data_ptr<float>()),
            cols4, rows);
    } else {
        dim3 grid((cols + block_size - 1) / block_size, grid_y);
        fused_ops_scalar_kernel<<<grid, block_size>>>(
            x_c.data_ptr<float>(),
            add_c.data_ptr<float>(),
            out.data_ptr<float>(),
            cols, rows);
    }

    return out;
}
"""

fused_ops_cpp_source = (
    "torch::Tensor fused_ops_cuda(torch::Tensor x, torch::Tensor add_val);"
)

# Compile the inline CUDA extension
fused_ops = load_inline(
    name="fused_activation_ops",
    cpp_sources=fused_ops_cpp_source,
    cuda_sources=fused_ops_source,
    functions=["fused_ops_cuda"],
    verbose=False,
    extra_cuda_cflags=["-O3"],
)


class ModelNew(nn.Module):
    """
    Optimized model: cuBLAS GEMM (addmm fuses the linear bias) followed by a
    single fused CUDA kernel for add_value + Swish + Tanh + GELU + Hardtanh.
    """
    def __init__(self, in_features, out_features, add_value_shape):
        super(ModelNew, self).__init__()
        self.matmul = nn.Linear(in_features, out_features)
        self.add_value = nn.Parameter(torch.randn(add_value_shape))

    def forward(self, x):
        # GEMM with fused bias epilogue via cuBLAS
        out = torch.addmm(self.matmul.bias, x, self.matmul.weight.t())
        # Single fused kernel for all remaining element-wise operations
        return fused_ops.fused_ops_cuda(out, self.add_value)
```